CN116817932B - Integrated method for track maintenance, track determination and real-time mapping of gravitational field - Google Patents
Integrated method for track maintenance, track determination and real-time mapping of gravitational field Download PDFInfo
- Publication number
- CN116817932B CN116817932B CN202211602711.7A CN202211602711A CN116817932B CN 116817932 B CN116817932 B CN 116817932B CN 202211602711 A CN202211602711 A CN 202211602711A CN 116817932 B CN116817932 B CN 116817932B
- Authority
- CN
- China
- Prior art keywords
- spacecraft
- thrust
- error
- coefficient
- law
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 80
- 238000013507 mapping Methods 0.000 title claims abstract description 41
- 238000012423 maintenance Methods 0.000 title claims abstract description 37
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 18
- 230000001133 acceleration Effects 0.000 claims description 71
- 239000011159 matrix material Substances 0.000 claims description 36
- 230000008859 change Effects 0.000 claims description 31
- 230000008569 process Effects 0.000 claims description 16
- 230000005484 gravity Effects 0.000 claims description 15
- 238000001914 filtration Methods 0.000 claims description 14
- 230000003044 adaptive effect Effects 0.000 claims description 11
- 238000011156 evaluation Methods 0.000 claims description 9
- 238000013461 design Methods 0.000 claims description 7
- 230000009471 action Effects 0.000 claims description 6
- 238000005259 measurement Methods 0.000 claims description 6
- 230000007704 transition Effects 0.000 claims description 6
- 230000033228 biological regulation Effects 0.000 claims description 4
- 230000000694 effects Effects 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 3
- 238000012360 testing method Methods 0.000 claims description 3
- 230000002708 enhancing effect Effects 0.000 abstract 1
- 238000004364 calculation method Methods 0.000 description 6
- 238000012546 transfer Methods 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000037406 food intake Effects 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Classifications
-
- 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
Landscapes
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
Abstract
The invention provides an integrated method for orbit maintenance, orbit determination and real-time mapping of a gravitational field, and belongs to the technical field of spacecraft orbit maintenance, orbit determination and star gravitational field mapping. The method comprises the following steps: s1, establishing a control law of track maintenance, a track determining algorithm and an estimation law of real-time mapping of a gravitational field; s2, tracking the set orbit of the spacecraft by the control law, an orbit determination algorithm and an estimation law. The invention solves the problems of deviation of the spacecraft from the expected orbit caused by errors of the thrust of the spacecraft, other unknown external perturbation forces and other factors, and has the advantages of improving the navigation precision, enhancing the control adaptability and robustness.
Description
Technical Field
The invention relates to the technical field of orbit maintenance and determination of a spacecraft and mapping of a gravitational field of a star, in particular to an integrated method for orbit maintenance, orbit determination and real-time mapping of the gravitational field.
Background
The trajectory of a spacecraft flying in space is called orbit. To accomplish a particular mission, the spacecraft needs to travel on a particular orbit. However, spacecraft are subject to space complications and tend to deviate from the desired orbit. For spacecrafts around stars, especially in deep space exploration, the factors affecting the orbit are mainly the high-order gravitational field of the stars, errors in orbit determination, uncertainty in the quality of the spacecraft itself, thrust errors of the spacecraft, and other unknown external ingestion forces. In order to maintain the spacecraft on the desired orbit, the spacecraft engine needs to be started to generate thrust, and the spacecraft is controlled to return to the desired orbit. This requires consideration of the three aspects of spacecraft orbit maintenance, orbit determination and star higher order gravitational field and issues.
For orbit maintenance problems, it is first necessary to study the orbit dynamics of a spacecraft in the gravitational field of stars. This is because the thrust of the existing spacecraft is very small compared with the gravity of the star, and the deviating effect of gravity on the orbit cannot be counteracted. Therefore, the existing methods are all used for researching the evolution law of spacecraft orbits under a star gravity field, and the orbits which are not easy to deviate and meet performance indexes, such as the frozen orbit of an earth satellite, the sun synchronous orbit and the like, are skillfully designed by utilizing an orbit evolution mechanism. Secondly, the orbit maintenance strategy of the spacecraft after being deviated from the orbit needs to be studied. The prior researches mostly use a track tracking method (a nominal track method), namely, a nominal track is designed in advance, then the real-time state (position, speed and the like) of the spacecraft is compared with the nominal track, and finally, the control law of the spacecraft for tracking the nominal track is designed. However, in the control law of the design orbit maintenance, the orbit deviation caused by the satellite high-order gravitational field is not considered in most cases.
For orbit determination, obtaining more accurate gravitational field data can improve navigation accuracy. The spacecraft near the earth can determine its own orbit through navigation constellations such as GPS, beidou, etc. Most spacecrafts flying in deep space can only rely on own satellite-borne equipment for navigation. The deep space exploration task is long in flight distance and long in flight time, the facing space environment is complex, and the ground measurement and control cannot meet the requirements of high precision and real-time performance. In addition, in a complex and changeable space environment, the failure rate of spacecraft components is greatly increased, and high-delay ground measurement and control cannot handle emergency situations. The navigation filtering process requires knowledge of the exact spacecraft dynamics, which requires knowledge of the high-order gravitational field data of the stars in advance or in real time.
For the problem of high-order gravity field, accurate gravity field data obtained by measurement is a precondition of utilizing the high-order gravity field. The process of measuring the high-order gravitational field is called mapping. CHAMP of DLR, GRACE with NASA and DLR in the united states, GOCE satellites in the aeronautical office, and GRAIL satellites in the united states make excellent contributions in mapping the earth's high-order gravitational field. The spacecraft runs through the annual and lunar orbits, and massive orbit and gravity gradient data are collected. These data are transmitted back to the earth and are subjected to long-term computation by high-performance computers or even supercomputers to obtain data of the earth's gravitational field. The traditional gravitational field mapping methods are all used for collecting data in advance and calculating on a computer afterwards, and have large calculated amount and no real-time property. However, in the detection task of the strange star, the calculation force of the spaceborne computer of the spacecraft is limited, the navigation data precision is not ideal, and therefore, the gravity field of the strange star is difficult to calculate by counting a large amount of orbit and gravity gradient data; meanwhile, the detection task hopes to acquire estimated data of the spherical harmonic coefficient in real time so as to gradually improve the track determination precision and ensure the success of track maintenance.
The invention provides an integrated strategy for track maintenance, track maintenance and real-time mapping of a gravitational field. The invention provides a concept of mapping a gravitational field, namely mapping the gravitational field of a strange star in real time in the orbital flight process. The mapping data are also applied to track determination and maintenance, so that navigation accuracy can be improved, and control adaptability and robustness can be enhanced.
Disclosure of Invention
The invention solves the technical problems that: when the spacecraft flies around a strange star, the high-order gravity field of the star is inaccurate, errors exist in orbit determination, the self mass of the spacecraft is uncertain, and the spacecraft deviates from an expected orbit due to factors such as errors in thrust of the spacecraft, other unknown external ingestion forces and the like.
In order to solve the problems, the technical scheme of the invention is as follows:
an integrated method for track maintenance, track determination and real-time mapping of gravitational field, comprising the following steps:
s1, establishing a control law of track maintenance, a track determining algorithm and an estimation law of real-time mapping of a gravitational field;
s2, tracking the set orbit of the spacecraft by the control law, an orbit determination algorithm and an estimation law.
Further, step S1 includes the steps of:
s1-1, establishing an orbit dynamics model
Modeling a satellite gravitational field by adopting a spherical harmonic coefficient method, and establishing a dynamic model of the spacecraft in the satellite harmonic coefficient gravitational field, wherein modeling conditions are as follows: during real-time mapping, spin axis direction of the star is unchanged, and the spacecraft mass is ignored relative to the target celestial mass;
s1-2, establishing a control law for calculating thrust required by the spacecraft from the current position and the current speed to the expected position and the expected speed, and automatically correcting a spherical harmonic coefficient by the spacecraft to achieve an estimation law for real-time mapping of a gravitational field;
S1-3, establishing an orbit determination algorithm for continuously calculating the current position and the current speed of the spacecraft to determine the current orbit of the spacecraft.
Further, step S2 includes the steps of:
s1-2-1, modeling thrust errors;
s1-2-2, designing a sliding mode surface;
s1-2-3, designing an estimation law;
s1-2-4, designing a control law;
s1-2-5, performing stability demonstration on the control law;
step S1-2-3, step S1-2-4, step S1-2-5 are performed simultaneously.
Further, step S1-2-1 includes the following:
the thrust of the thruster for the spacecraft is provided: neglecting the inertia of the system, namely considering that the control action is applied without time delay, and assuming that the thrust of the thruster has no error in the direction of the force and only has error in the magnitude, establishing a model formula of the thrust error of the thruster as follows:
in the above formula, u is the actual thrust output by the thruster, u t For theoretical thrust calculated by control law, ε is the relative error of thrust, ε 1 As the first coefficient of the thruster, epsilon 2 For the second coefficient of the thruster, m is the mass of the spacecraft, and it is assumed that in ground tests the upper bound epsilon of the relative error epsilon of the thrusts can already be known max And is also provided with
Further, step S1-2-2 includes the following:
according to the perturbation force caused by the nominal model attraction force and nominal model error of the spacecraft in the celestial body gravitational field and the control force of the spacecraft self-propulsion system, a kinetic equation under the nominal model is established as follows:
In the above-mentioned method, the step of,is the acceleration of the spacecraft, f is the gravitational acceleration under a nominal model, D is the gravitational acceleration caused by model errors, and +.>u is the actual thrust output by the thruster, m is the mass of the spacecraft, and the definition of the gravitational acceleration f is as follows:
f＝T -1 (ω pr t)f b
in the above, f is gravity acceleration, f b Is the projection of f in the star fixedly connected coordinate system, omega pr Is an a priori estimate of the rotational angular velocity of the star,wherein omega pr Is the rotational angular velocity of the star, t is the time,
the position control error formula of the spacecraft is defined as follows:
e 1 ＝r-χ d
defining the speed control error of the spacecraft as follows:
in the above, e 1 E is the position control error of the spacecraft 2 Is the position control error of the spacecraft, r is the position of the spacecraft,for spacecraft speed χ d ∈R 3 For nominal position +.>At the time of the nominal speed of the vehicle,
the definition formula of the slip plane s is:
in the above, S is a slip form surface, S 1 Is a position error, S 2 For speed error, p is the first coefficient of sliding mode surface, the value is positive constant, q is the second coefficient of sliding mode surface, the value is normal number, and the control law is implementedThe target of (c) becomes s=0.
Further, step S1-2-3 includes the following:
the adaptive estimation is realized by adopting a method of correcting parameters in real time, the perturbation acceleration caused by a nominal dynamics model relative to a real model is bounded, and the upper bound of the perturbation acceleration is set as D max ∈R 3 ，D≤D max The mass m of the spacecraft and the first coefficient epsilon of the thruster are calculated 1 First coefficient epsilon of thruster 2 And upper boundary of perturbation acceleration D max For the adaptive estimator, an estimation error is defined as:
in the above-mentioned method, the step of,error of estimated value for spacecraft quality, +.>Is the estimated value of the quality of the spacecraft, m is the quality of the spacecraft,error of upper bound estimate for perturbation acceleration, +.>For the upper bound estimate of the perturbation acceleration, D max For the upper bound of perturbation acceleration, +.>Error of estimated value of spherical harmonic coefficient of gravitational field of star, < +.>Is the estimated value of the harmonic coefficient of the gravitational field of the star lambda j Is the gravitational field spherical harmonic coefficient of the star, lambda 1 ＝C 20 ，λ 2 ＝C 21 Sequentially numbered, the number is that,
the adaptive estimation law definition formula is:
in the above-mentioned method, the step of,the change rate of the error of the spacecraft quality estimated value is represented by q, the second coefficient of the sliding mode surface, gamma, the first estimated law parameter and m 0 ＝m 0 (t) is the upper mass limit of the spacecraft, and the condition that at any moment t, m is satisfied 0 (t)≥m(t)/>For spacecraft quality estimation, S is the sliding mode surface, sgn is the sign function, ++>Is the difference between the true acceleration and the nominal angular velocity, +.>For estimating the rate of change of the thruster coefficient, E i For the second estimated law parameter, Ω is the third estimated law parameter, < >>For the rate of change of the perturbation acceleration upper bound estimate error, Λ D For the fourth estimated law parameter, γ, Ω, E i 、Λ D All are normal numbers and are added with->For the rate of change of the spherical harmonic coefficient estimation value error, Λ λj For the fifth estimated law parameter D λj Is the spherical harmonic coefficient lambda j The resulting gravitational acceleration perturbation,
obtaining each spherical harmonic coefficient lambda j The estimated value and the estimated value of the star autorotation speed omega are obtained in real time, namely the star gravitational field and the Coriolis acceleration caused by the star autorotation are obtained in real time, namely the real-time mapping of the star gravitational field is realized,
D λj about lambda j Partial derivative d of (2) λj The solution formula of (2) is:
this can be achieved by:
in the above-mentioned method, the step of,for the rate of change of the spherical harmonic coefficient estimation value error, Λ λj For the fifth estimated law parameter d λj For D λj About lambda j Q is the second coefficient of the sliding mode surface, m 0 For the upper mass boundary of the spacecraft, S is the sliding mode surface, sgn is the sign function, ++>For the spherical harmonic coefficient estimation value, < >>To estimate the thrust, u t For theoretical thrust->For the estimation of the thruster coefficients, < >>For the first thruster coefficient estimation, +.>For the second thruster coefficient estimation, +.>As the estimated value of the quality of the spacecraft,
due to the actual thrust u and the estimated thrustIs smaller and the thruster is specific to the pulse I sp Very large, it follows that:
in the above-mentioned method, the step of,g is the real mass change rate of the spacecraft 0 For the gravitational acceleration of the earth's surface, I sp For thrust specific thrust, u is the actual thrust, +. >In order to estimate the thrust force,
the rate of change of the spacecraft mass estimate can be obtained as:
in the above-mentioned method, the step of,for the rate of change of the quality estimate +.>Error rate of change for spacecraft mass estimation, < >>For the rate of change of the real mass, m 0 Is the upper mass limit of the spacecraft, S is a sliding mode surface, sgn is a sign function, g 0 For the gravitational acceleration of the earth's surface, I sp Is thrust device specific thrust->To estimate the thrust +.>Q is a second coefficient of the sliding mode surface, and gamma is a first estimation law parameter.
Still further, step S1-2-4 includes the following:
the control law expression is as follows:
in the above, u t As theoretical thrust, u t1 As the first term of theoretical thrust, u t2 S is a sliding mode surface, S is a theoretical thrust second term 1 Is a position error, S 2 P is a first coefficient of a sliding mode surface and is a positive constant; q is the second coefficient of the sliding die surface, and m is the normal number 0 For the upper mass limit of the spacecraft, k D For a first thrust adjustment factor, k D Not less than 1, sgn is a sign function,for the upper bound estimate of the perturbation acceleration, k ε For a second thrust adjustment factor, k ω For a third thrust adjustment factor, k s For the fourth thrust adjustment factor,for the estimated value of the first thruster coefficient, +. >For the estimated value of the second thruster coefficient, +.>Is the estimated value of the rotation speed of the star body, D λj Is the spherical harmonic coefficient lambda j Induced gravitational acceleration perturbation, ε max Is the upper limit value of the relative error epsilon of the thrust forcek ε 、k ω And k s Is positive constant for adjusting epsilon i Omega and lambda j Impact on thrust, greater k ε 、k ω And k s Will decrease epsilon i Omega and lambda j Too little effect, undermining the control law's adaptivity,
the formula for defining the sign function sgn is:
sgnv＝diag[sgn(v 1 ) sgn(v 2 ) sgn(v 3 )]
in the above formula, v is any vector, v 1 V is the first component of v 2 V is the second component of v 3 Is the third component of v.
Preferably, step S1-2-5 comprises the following:
the function for stability demonstration was constructed as follows:
V＝V 0 +V 1
in the above formula, V is a stability proving function, V 0 For the first proof of stability function, V 1 For the second proof function of stability, m is the spacecraft mass, S is the sliding mode surface, gamma is the first estimation law parameter,k is the estimated quantity of the quality of the spacecraft ε For a second thrust adjustment factor, k ω For a third thrust adjustment factor, k s For a fourth thrust regulation factor, E i For the second evaluation law parameter,/>For the i-th estimated value of the thruster coefficient, Ω is the third estimated law parameter, ++>Is the estimated value of the rotation angular velocity of the star, Λ λj For the fifth evaluation law parameter, ++ >Is the spherical harmonic coefficient lambda j Induced gravitational acceleration perturbation, Λ D For the third evaluation law parameter, ++>For the error of the upper bound estimate of the perturbation acceleration,
as can be seen, all the terms of the stability proving function V are smaller than or equal to 0, soAnd->The filling condition of (2) is s=0, and the LaSalle invariant set theorem can know that the system converges global asymptote to s=0, namely s 1 =0 and s 2 =0, both the position error and the velocity error converge to zero.
Preferably, step S1-3 includes the following:
under the action of the gravitational field of the star and the thrust of the star, the continuous state equation of the spacecraft is as follows:
in the above formula, X (t) is a state vector and X (t) = [ X y z v ] x v y v z ] T U is the thrust of the thruster and u= [ u ] x u y u z ]M is the mass of the spacecraft, W (t) is the error of the state equation,
linearizing the continuous state equation by setting a sampling step length to obtain a discretized continuous state equation:
X k ＝F k X k-1 +B k u k +w k
in the above, X k X is the state of spacecraft at time k k-1 Is the state of the spacecraft at the moment k-1, F k For state transition matrix, B k For the control matrix at time k, u k To control vector, w k Is state noise and wk= [ w k，x w k，y w k，z w k，vx w k，vy w k，vz ]From errors caused by inaccurate estimated values of spherical harmonic coefficients, spacecraft mass and thruster coefficients,
the observation equation for spacecraft position is:
In the above, Z k Is the observed value of the current position of the spacecraft, v k To observe noise H k In order to observe the matrix,
the spacecraft is set to directly observe the position of the spacecraft to obtain an observation matrix H of the spacecraft k The following are provided:
assumed state noise w k And observation noise v k Is zero-mean Gaussian noise, independent of each other, and let Q k ＝cov(w k )，R k ＝cov(v k ) The recurrence equation for the extended kalman filter is obtained as:
in the above-mentioned method, the step of,to estimate the state at time k-1 for time k, F k To act on the state transition matrix at time k-1, X k-1，k-1 To estimate the state at time k-1, u k To control vector B k To act on the control vector u k Control matrix, w k P is process noise k，k-1 For the state prediction error covariance matrix from the moment k-1 to the moment k, P k-1，k-1 As the estimated error covariance matrix at the moment k-1, Q k K is a process noise covariance matrix k For Kalman gain, H k For the observation equation, R k To observe covariance matrix of noise, X k，k As an estimated value of the state at time k, Z k P is the measurement value of the state at the moment k k，k As an estimation error covariance matrix at the k moment,
setting the prior value error to be-20%, and obtaining the spherical harmonic coefficient input value of the extended Kalman filter as follows:
in the above, lambda i，put For the spherical harmonic coefficient values for kalman filtering,lambda is the estimated value of the spherical harmonic coefficient p Is the prior value of the spherical harmonic coefficient.
Further preferably, step S2 comprises the steps of:
s2-1, determining a set orbit of the spacecraft operation by a navigation task of the spacecraft, wherein the set orbit comprises expected positions and expected speeds of the spacecraft at all moments;
s2-2, continuously calculating the current position and the current speed of the spacecraft through an orbit determination algorithm in the navigation process of the spacecraft, and determining the current expected position and the expected speed through a set orbit;
s2-3, when errors exist in the current position and the current speed of the spacecraft, the expected position and the expected speed of the set orbit, the control law calculates the thrust of the spacecraft internal thruster required by the spacecraft from the current position and the current speed to the expected position and the expected speed, and in the control law calculation process, the estimation law provides the estimation value of the parameter required by the control law;
s2-4, the spacecraft controls the thruster to achieve the thrust calculated by the control law in the step S2-3, so that the spacecraft achieves the expected position and the expected speed of the set orbit, and the tracking of the spacecraft on the set orbit is achieved.
The beneficial effects of the invention are as follows:
(1) The mapping method provided by the invention is real-time, has less calculation time and less calculation amount compared with the traditional mapping method, and can be applied to a filter for track determination, so that the navigation precision can be improved, the method is applied to a control law for track maintenance, and the control adaptability and robustness can be enhanced;
(2) The invention provides a control design method for changing the thrust, which introduces coefficients into a sliding mode surface, adjusts the thrust by adjusting the sliding mode surface coefficients in the control process, and can reduce the thrust by reducing the sliding mode surface coefficients when the control early-stage error is very large, so as to avoid exceeding the maximum thrust; in the middle control stage, the thrust can be increased by increasing the sliding mode surface coefficient, so that the error convergence speed is increased.
Drawings
FIG. 1 is a flow chart of an integrated method of orbit maintenance, orbit determination and real-time mapping of gravitational fields in accordance with an embodiment;
FIG. 2 is a flowchart of the method steps S1-2 of integrating rail maintenance, rail determination and real-time mapping of gravitational field in accordance with an embodiment;
FIG. 3 is a diagram depicting the position of a spacecraft in two coordinate systems in an embodiment;
FIG. 4 is a schematic diagram of extended Kalman filtering in an embodiment;
FIG. 5 is a logical relationship of orbit maintenance, gravitational field mapping, autonomous navigation in an embodiment.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail below with reference to the accompanying drawings, and it is apparent that the described embodiments are only some embodiments of the present invention, not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
The terminology used in the embodiments of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used in this application and the appended claims, the singular forms "a," "an," and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise, the "plurality" generally includes at least two.
It should be understood that although the terms first, second, third, etc. may be used to describe … … in embodiments of the present invention, these … … should not be limited to these terms. These terms are only used to distinguish … …. For example, the first … … may also be referred to as the second … …, and similarly the second … … may also be referred to as the first … …, without departing from the scope of embodiments of the present invention.
Example 1
The embodiment is an integrated method for track maintenance, track determination and real-time mapping of a gravitational field, comprising the following steps:
s1, establishing a control law of track maintenance, a track determining algorithm and an estimation law of real-time mapping of a gravitational field, wherein the method comprises the following steps of:
s1-1, establishing an orbit dynamics model
Assuming that the spin axis orientation of the star is unchanged during real-time mapping, i.e., the star rotates about its own maximum inertia axis, the spacecraft mass is negligible with respect to the target celestial mass,
Modeling a satellite gravitational field by adopting a spherical harmonic coefficient method, and establishing a dynamic model of the spacecraft in the satellite harmonic coefficient gravitational field, wherein the method comprises the following steps of:
s1-1-1, and assumption conditions for establishing orbit dynamics model
Because the mission time of real-time mapping is much smaller than the change time of spin axis direction of the star, it is assumed that the spin axis direction of the star is unchanged during real-time mapping, i.e. the star rotates around its own maximum inertia axis, and at the same time, it is assumed that the spacecraft mass is negligible with respect to the target celestial mass,
s1-2, and establishing a coordinate system
Establishing a fixed coordinate system O-XYZ and an inertial coordinate system O-XYZ, wherein O is the mass center of the strange star, the Z axis and the X axis of the fixed coordinate system are respectively along the maximum rotational inertia axis and the minimum rotational inertia axis of the star, the Y axis is determined by a right hand rule, the Z axis of the inertial coordinate system is along the maximum inertia main axis of the star, the X axis and the Y axis are respectively overlapped with the X axis and the Y axis at the initial moment,
the positions of the spacecraft in the fixed coordinate system O-XYZ and the inertial coordinate system O-XYZ are respectively expressed as R= [ X Y Z ]]T and r= [ x y z ]] T As shown in fig. 3, the present invention,
the transfer matrix from the inertial coordinate system to the fixed coordinate system is:
in the above formula, T is a transfer matrix from an inertial coordinate system to a fixedly connected coordinate system, ω is the star rotation angular velocity, T is time,
Three spherical coordinates of the spacecraft under the fixed coordinate system are defined as follows: r represents the distance from the spacecraft to the origin O, phi represents latitude, eta represents longitude, and the conversion formula of the fixed coordinate system and the inertial coordinate is as follows:
R＝||R|||.
in the above, X is the X-axis component of the coordinate of the spacecraft in the fixed coordinate system O-XYZ, Y is the Y-axis component of the coordinate of the spacecraft in the fixed coordinate system O-XYZ, Z is the Z-axis component of the coordinate of the spacecraft in the fixed coordinate system O-XYZ,
s1-1-3, and establishing an orbit dynamics model
Modeling the gravitational field of the star by using a spherical harmonic coefficient method, wherein the calculation formula of the gravitational potential of the check point is as follows:
in the above formula, U (R) is the electromotive force of the inspection point, R e For reference radius of star, C nm Is the first spherical harmonic coefficient, S nm Is the second spherical harmonic coefficient, P nm (u) is a Legendre function, P n (u) is a Legendre polynomial,
first term spherical harmonic coefficient C nm The expression formula of (2) is as follows:
second spherical harmonic coefficient S nm The expression formula of (2) is as follows:
legendre function P nm The expression formula of (u) is as follows:
legendre polynomial P n The expression formula of (u) is as follows:
in the above, P n (u) is a Legendre polynomial,
the known value of the spherical harmonic coefficient, i.e. the first spherical harmonic coefficient C nm Second spherical harmonic coefficient S nm Legendre function P nm (u), legendre polynomial P n (u) substitution into the force potential, i.eNominal model U from which the induced potential can be derived 0 Nominal model of the induced force U 0 The relation with the true attraction U is:
U＝U 0 +U D
in the above, U D For nominal model U 0 The difference from the true gravitational potential U,
the gravitational acceleration of the spacecraft in the gravitational field of the star natural model is as follows:
in the above, a is gravitational acceleration of a spacecraft in a gravitational field of a star natural model, U is true gravitational potential, U 0 As a nominal model of the induced force, U D For nominal model U 0 The difference from the true gravitational potential U, x is the x-axis component of the position of the spacecraft in the inertial coordinate system O-xyz, y is the y-axis component of the position of the spacecraft in the inertial coordinate system O-xyz, and z is the z-axis component of the position of the spacecraft in the inertial coordinate system O-xyz;
s1-2, establishing a control law for calculating thrust required by the spacecraft from the current position, the current speed to the expected position and the expected speed, and automatically correcting the spherical harmonic coefficient by the spacecraft to achieve an estimation law for real-time mapping of a gravitational field, wherein the method comprises the following steps of:
s1-2-1, thrust error modeling
The thrust of the thruster for the spacecraft is provided: neglecting the inertia of the system, namely considering that the control action is applied without time delay, and assuming that the thrust of the thruster has no error in the direction of the force and only has error in the magnitude, establishing a model formula of the thrust error of the thruster as follows:
In the above formula, u is the actual thrust output by the thruster, u t For theoretical thrust calculated by control law, ε is the relative error of thrust, ε 1 Is a thrusterFirst coefficient, epsilon 2 For the second coefficient of the thruster, m is the mass of the spacecraft, and it is assumed that in ground tests the upper bound epsilon of the relative error epsilon of the thrusts can already be known max And is also provided with
S1-2-2, slip form surface design
According to the nominal model gravitation of the spacecraft in the celestial body gravitational field, the perturbation force caused by the model error and the control force of the spacecraft self-propulsion system, a kinetic equation under the nominal model is established as follows:
in the above-mentioned method, the step of,the method is characterized in that the method comprises the steps of (1) the acceleration of a spacecraft, f is the gravitational acceleration under a nominal model, D is the gravitational acceleration caused by model errors, u is the actual thrust output by a thruster, m is the mass of the spacecraft, and the gravitational acceleration f is defined as follows:
f＝T -1 (ω pr t)f b
wherein f b Is the projection of f in the star fixedly connected coordinate system, omega pr Is an a priori estimate of the rotational angular velocity of the star,wherein omega pr The rotation angular velocity of the star is given, and t is time.
The position control error formula of the spacecraft is defined as follows:
e 1 ＝r-χ d
defining the speed control error of the spacecraft as follows:
in the above, e 1 E is the position control error of the spacecraft 2 Is the position control error of the spacecraft, r is the position of the spacecraft, For spacecraft speed χ d ∈R 3 For nominal position +.>At the time of the nominal speed of the vehicle,
the aim of the control is therefore to maintain the spacecraft stable on the nominal orbit, i.e. to control the position and velocity errors to 0, i.e. e 1 ＝0，e 2 =0, substituting the control error definition into the kinetic equation yields:
wherein,
the definition formula of the slip plane s is:
in the above, S is a slip form surface, S 1 Is a position error, S 2 For the speed error, p is the first coefficient of the sliding mode surface, the value is a positive constant, q is the second coefficient of the sliding mode surface, the value is a positive constant, p and q are real numbers larger than zero, and s is because 1 、s 2 Is a time-varying quantity, p and q are variable constants, and it is considered that the sufficient requirement for s=0 is s 1 =0 and s 2 =0. The target of the control law becomes s=0.
S1-2-3, design of estimation law
The adaptive control can automatically adjust the characteristics of the system itself to solve the problems of the variation of the uncertain object and the dynamic characteristics thereof, and is mainly classified into an indirect type and a direct type. The indirect self-adaptive control requires real-time identification of system parameters, and proper control is designed based on the real-time identification, and the direct self-adaptive control does not need to identify the system parameters but directly updates the control system parameters, and the embodiment adopts a method of correcting the parameters in real time to realize self-adaptive estimation:
The adaptive estimation is realized by adopting a method of correcting parameters in real time, the perturbation acceleration is bounded due to the error of a nominal dynamics model relative to a real model, and the upper bound of the perturbation acceleration is set as D max ∈R 3 Is a constant vector with positive values of all components, and meets the following conditions:
D x ≤D max，x D y ≤D max，y D z ≤D max，z
simply denoted as D is less than or equal to D max First coefficient epsilon of thruster 1 First coefficient epsilon of thruster 2 And upper boundary of perturbation acceleration D max For the adaptive estimator, an estimation error function is defined as:
in the above-mentioned method, the step of,error of estimated value for spacecraft quality, +.>Is the estimated value of the quality of the spacecraft, m is the quality of the spacecraft,error of upper bound estimate for perturbation acceleration, +.>For the upper bound estimate of the perturbation acceleration, D max For the upper bound of perturbation acceleration, +.>Error of estimated value of spherical harmonic coefficient of gravitational field of star, < +.>Is the estimated value of the harmonic coefficient of the gravitational field of the star lambda j Is the gravitational field spherical harmonic coefficient of the star, lambda 1 ＝C 20 ，λ 2 ＝C 21 Sequentially numbered, the number is that,
the adaptive estimation law definition formula is:
/>
in the above，The change rate of the error of the spacecraft quality estimated value is represented by q, the second coefficient of the sliding mode surface, gamma, the first estimated law parameter and m 0 ＝m 0 (t) is the upper mass limit of the spacecraft, and the condition that at any moment t, m is satisfied 0 (t)≥m(t)/>For spacecraft quality estimation, S is the sliding mode surface, sgn is the sign function, ++ >Is the difference between the true acceleration and the nominal angular velocity, +.>For estimating the rate of change of the thruster coefficient, E i For the second estimated law parameter, Ω is the third estimated law parameter, < >>For the rate of change of the perturbation acceleration upper bound estimate error, Λ D For the fourth estimated law parameter, γ, Ω, E i 、Λ D All are normal numbers and are added with->For the rate of change of the spherical harmonic coefficient estimation value error, Λ λj For the fifth estimated law parameter D λj Is the spherical harmonic coefficient lambda j The resulting gravitational acceleration perturbation,
obtaining each spherical harmonic coefficient lambda j The estimated value and the estimated value of the star autorotation speed omega are obtained in real time, namely the star gravitational field and the Coriolis acceleration caused by the star autorotation are obtained in real time, namely the real-time mapping of the star gravitational field is realized,
D λj about lambda j Partial derivative d of (2) λj The solution formula of (2) is:
this can be achieved by:
in the above-mentioned method, the step of,for the rate of change of the spherical harmonic coefficient estimation value error, Λ λj For the fifth estimated law parameter d λj For D λj About lambda j Q is the second coefficient of the sliding mode surface, m 0 For the upper mass boundary of the spacecraft, s is the slip plane, sgn is the sign function, ++>For the spherical harmonic coefficient estimation value, < >>To estimate the thrust, u t For theoretical thrust->For the estimation of the thruster coefficients, < >>For the first thruster coefficient estimation, +.>For the second thruster coefficient estimation, +. >As the estimated value of the quality of the spacecraft,
due to the actual thrust u and the estimated thrustIs smaller and the thruster is specific to the pulse I sp Very large, it follows that:
in the above-mentioned method, the step of,g is the real mass change rate of the spacecraft 0 For the gravitational acceleration of the earth's surface, I sp For thrust specific thrust, u is the actual thrust, +.>In order to estimate the thrust force,
the rate of change of the spacecraft mass estimate can be obtained as:
/>
in the above-mentioned method, the step of,for the rate of change of the quality estimate +.>Error rate of change for spacecraft mass estimation, < >>For the rate of change of the real mass, m 0 Is the upper mass limit of the spacecraft, S is a sliding mode surface, sgn is a sign function, g 0 For the gravitational acceleration of the earth's surface, I sp Is thrust device specific thrust->To estimate the thrust +.>Q is the second coefficient of the sliding mode surface, gamma is the first estimation law parameter,
note that in the estimation law, the gravitational field sets two values: the spherical harmonic coefficient and the upper bound of the perturbation acceleration are theoretically enough, and the upper bound of the perturbation acceleration is not required to be estimated in the recovery process. In practice we can use the spherical harmonic coefficients and upper bound of the perturbation acceleration to adjust the thrust for better control. However, if some estimates of the spherical harmonic coefficients and the star self-rotational angular velocity are miscalculated, the thrust will be biased and even out of range. This may lead to failure of track maintenance. Therefore, it is necessary to estimate the upper bound of the perturbation acceleration and introduce it into the thrust. In the presence of control errors, the estimate will increase until track maintenance is achieved. Thus, by introducing an estimate of the upper bound of the perturbation acceleration, the robustness of the control is improved. At the end of the flight, the estimated values of these parameters will be close to the true values. But in the later stages, navigation errors are constantly present, so that the control errors cannot converge to zero. This may result in the estimated value changing continuously, rather than being stable at the true value. Therefore, a stop condition needs to be set to the estimation process. The stopping criterion set by this patent is a position tracking error, and if the error falls below a certain threshold and is maintained for a period of time, the estimation is considered to be stopped. After the estimated value is obtained, the spacecraft can design a new fuel-saving orbit for the following tasks by utilizing dynamics near the new problem. For a transfer from a previous track to a new track, it is only necessary to define the new track as a nominal track.
S1-2-4, control law design
The control law expression is as follows:
in the above, u t As theoretical thrust, u t1 As the first term of theoretical thrust, u t2 S is a sliding mode surface, S is a theoretical thrust second term 1 Is a position error, S 2 P is a first coefficient of a sliding mode surface and is a positive constant; q is the second coefficient of the sliding die surface, and m is the normal number 0 For the upper mass limit of the spacecraft, k D For a first thrust adjustment factor, k D Not less than 1, sgn is a sign function,for the upper bound estimate of the perturbation acceleration, k ε For a second thrust adjustment factor, k ω For a third thrust adjustment factor, k s For the fourth thrust adjustment factor,for the estimated value of the first thruster coefficient, +.>For the estimated value of the second thruster coefficient, +.>Is the estimated value of the rotation speed of the star body, D λj Is the spherical harmonic coefficient lambda j Induced gravitational acceleration perturbation, ε max Is the upper limit value of the relative error epsilon of the thrust forcek ε 、k ω And k s Is positive constant for adjusting epsilon i Omega and lambda j Impact on thrust, greater k ε 、k ω And k s Will decrease epsilon i Omega and lambda j Too little effect, undermining the control law's adaptivity,
the formula for defining the sign function sgn is:
sgnv＝diag[sgn(ν 1 ) sgn(v 2 ) sgn(v 3 )]
in the above formula, v is any vector, v 1 A first component of v, a second component of v2, a third component of v3,
In order to improve the control efficiency and reduce the thrust as much as possible, the principle of selecting the nominal state on the nominal track is as follows: on the nominal track, the state of the position closest to the current real position is the nominal state.
Will s=ps 1 +qs 2 Substitution into the control law can be:
it can be seen that the p/q value determines the position error s 1 And speed error s 2 The specific gravity of the thrust force. In this equation, if the q value is unchanged and p is smaller, the position error s 1 The smaller the specific gravity in the thrust force, the speed error s 2 The greater the specific gravity. Because the estimation law contains the parameter q, the thrust is adjusted by adopting a method of adjusting p by fixing q for the convenience of parameter adjustment in the estimation law. In the early control stage, the spacecraft has larger position error, which possibly causes overlarge thrust and exceeds the maximum thrust which can be provided by the thruster, and the thrust can be reduced by reducing the p value; in the middle of the control, the position error is small, and the thrust force can be increased by increasing the p value.
In this way, by adjusting the p value in the middle and front stages of control, the saturation thrust can be maintained, and the convergence speed can be increased. To maintain the saturation thrust, a suitable p value is obtained, and a dichotomy may be used. In the latter stage of control, the speed error s can be increased by decreasing the p value 2 In thrust forceThe proportion of the track is controlled by differential control, thereby tracking the nominal track more stably. When the system moves near the sliding mode surface, repeated sign changes of the sign function sgn occur, and high-frequency flutter is caused. It is the control force that is subject to high frequency magnitude abrupt changes. The control device cannot handle such high frequency chatter, and may be damaged, causing control failure. In order to be able to attenuate the flutter of the thrust, the sign function sgnv is replaced by the satv function:
satv＝diag[sat(ν 1 ) sat(v 2 ) sat(v 3 )]
wherein delta is i I.e. the boundary layer thickness, can be set as a constant or as an adaptive variable. The control accuracy and the boundary layer thickness delta of the method i In this regard, the smaller the thickness, the higher the accuracy.
Another approach is to replace sgnv with a smoother hyperbolic tangent tanh v function, namely:
tanhv＝diag[tanh(a 1 ν 1 ) tanh(a 2 v 2 ) tanh(a 3 v 3 )]
wherein a is 1 、a 2 、a 3 All are normal numbers.
S1-2-5, and stability of control law is proved
The function for stability demonstration was constructed as follows:
V＝V 0 +V 1
in the above formula, V is a stability proving function, V 0 For the first proof of stability function, V 1 For the second proof function of stability, m is the spacecraft mass, S is the sliding mode surface, gamma is the first estimation law parameter,k is the estimated quantity of the quality of the spacecraft ε For a second thrust adjustment factor, k ω For a third thrust adjustment factor, k s For a fourth thrust regulation factor, E i For the second evaluation law parameter,/>For the i-th estimated value of the thruster coefficient, Ω is the third estimated law parameter, ++>Is the estimated value of the rotation angular velocity of the star, Λ λj For the fifth evaluation law parameter, ++>Is the spherical harmonic coefficient lambda j Induced gravitational acceleration perturbation, Λ D For the third evaluation law parameter, ++>For the error of the upper bound estimate of the perturbation acceleration,
V 0 deriving time t and taking the s derivative of t as:
on the other hand, V 1 For time of dayt is derived to obtain:
the derivative of the function V with respect to time t is therefore:
obviously, can obtainThe first and second terms of (2) are equal to or less than 0:
and because of m 0 ＞m，D＜D max Constant is true, so the following is obtained:
qs T [mD-m 0 sgns·D max ]≤0
further, the method comprises the steps of,hengzhi and->With upper bound, let k D 、k s Take a sufficiently largeThe value is obtained by:
as can be seen, all the terms of the stability proving function V are smaller than or equal to 0, soAnd->The filling condition of (2) is s=0, and the LaSalle invariant set theorem can know that the system converges global asymptote to s=0, namely s 1 =0 and s 2 =0, both the position error and the velocity error converge to zero,
step S1-2-3, step S1-2-4 and step S1-2-5 are carried out simultaneously;
s3, establishing an orbit determination algorithm for continuously calculating the current position and the current speed of the spacecraft to determine the current orbit of the spacecraft, wherein the orbit determination algorithm comprises the following contents:
The orbit determination method of the embodiment comprehensively considers the astronomical observation value and orbit dynamics, and obtains the optimal estimation of navigation information through filtering. In the filtering method, extended Kalman filtering is adopted, fig. 4 shows a process of orbit determination, an observation device such as a star sensor and an optical camera obtains an original observation value, an observation model is input to obtain a navigation observation value, the spacecraft quality, thrust and spherical harmonic coefficient estimation value obtained by a parameter estimation module is input to a state model, a navigation prediction value based on orbit dynamics is obtained, and the observation value and the prediction value based on dynamics are subjected to extended Kalman filtering and output as an optimal estimation value of the navigation state of the spacecraft.
Under the action of the gravitational field of the star and the thrust of the star, the continuous state equation of the spacecraft is as follows:
in the above formula, X (t) is a state vectorQuantity and X (t) = [ X y z v ] x v y v z ] T U is the thrust of the thruster and u= [ u ] x u y u z ]M is the mass of the spacecraft, W (t) is the error of the state equation, in practical application,an estimate of the spherical harmonic coefficients needs to be used, therefore +.>An error exists between the estimated value and the actual value of (a); />The middle thrust and the mass are also the estimated thrust used +.>And estimate quality->There is also an error between its value and the actual value. The error W (t) of the state equation is therefore used to characterize the error due to the spherical harmonic coefficients and the inaccuracy of the spacecraft mass, thruster coefficient estimates.
Linearizing the continuous state equation by setting a sampling step length to obtain a discretized continuous state equation:
X k ＝F k X k-1 +B k u k +w k
in the above, X k X is the state of spacecraft at time k k-1 Is the state of the spacecraft at the moment k-1, F k For state transition matrix, B k For the control matrix at time k, u k To control vector, w k Is state noise and w k ＝[w k，x w k，y w k，z w k，vx w k，vy w k，vz ] T From errors caused by inaccurate estimated values of spherical harmonic coefficients, spacecraft mass and thruster coefficients,
the observation equation for spacecraft position is:
in the above, Z k Is the observed value of the current position of the spacecraft, v k To observe noise H k In order to observe the matrix,
the emphasis of the autonomous navigation method is that estimated values of parameters such as spacecraft quality, thrust, star harmonic coefficients and the like are introduced into an autonomous navigation algorithm, so that the error of a navigation observation value is corrected. What kind of observation means is not the key point of the article, so the simplest observation means is selected, and the spacecraft is set to directly observe the position of the spacecraft to obtain the observation matrix H of the spacecraft k The following are provided:
assumed state noise w k And observation noise v k Is zero-mean Gaussian noise, independent of each other, and let Q k ＝cov(w k )，R k ＝cov(v k ) As the control proceeds, each estimated value becomes accurate gradually, so the state noise w k Should decrease gradually as the estimation proceeds.
The recurrence equation for the extended kalman filter is obtained as:
in the above-mentioned method, the step of,to estimate the state at time k-1 for time k, F k To act on the state transition matrix at time k-1, X k-1，k-1 To estimate the state at time k-1, u k To control vector B k To act on the control vector u k Control matrix, w k P is process noise k，k-1 For the state prediction error covariance matrix from the moment k-1 to the moment k, P k-1，k-1 As the estimated error covariance matrix at the moment k-1, Q k For process noise covariance matrix, K k For Kalman gain, H k For the observation equation, R k To observe covariance matrix of noise, X k，k As an estimated value of the state at time k, Z k P is the measurement value of the state at the moment k k，k As an estimation error covariance matrix at the k moment,
F k the matrix contains gravitational accelerationRelative to r k-1 I.e. the spherical harmonic coefficient information of the gravitational field of the star must be utilized. In the simulation process, the priori value of the star fourth-order spherical harmonic coefficient is considered to exist, and certain errors exist, and the error range is known. Therefore, by comparing the prior value and the estimated value of the spherical harmonic coefficient, the value closest to the true value can be selected from the prior value and the estimated value, and the value is input into the EKF, so that the navigation accuracy is improved. Setting the prior value error to be-20%, and obtaining the spherical harmonic coefficient input value of the extended Kalman filter as follows:
In the above, lambda input For the spherical harmonic coefficient values for kalman filtering,lambda is the estimated value of the spherical harmonic coefficient p Is the prior value of the spherical harmonic coefficient.
S2, tracking the set orbit of the spacecraft by the control law, an orbit determination algorithm and an estimation law, wherein the method comprises the following steps of:
s2-1, determining a set orbit of the spacecraft operation by a navigation task of the spacecraft, wherein the set orbit comprises expected positions and expected speeds of the spacecraft at all moments,
s2-2, continuously calculating the current position and the current speed of the spacecraft through an orbit determination algorithm in the navigation process of the spacecraft, determining the current expected position and the expected speed through a set orbit,
s2-3, when errors exist between the current position and the current speed of the spacecraft and the expected position and the expected speed of the set orbit, the control law calculates the thrust of the thruster in the spacecraft required by the spacecraft from the current position and the current speed to the expected position and the expected speed, and in the control law calculation process, the estimation law provides the estimation value of the parameter required by the control law,
s2-4, the spacecraft controls the thruster to achieve the thrust calculated by the control law in the step S2-3, so that the spacecraft achieves the expected position and the expected speed of the set orbit, and the tracking of the spacecraft on the set orbit is achieved.
As shown in fig. 1, the estimation law estimates the detector mass, the thruster coefficient, the star self-rotation angle speed and the spherical harmonic coefficient according to the tracking error and the current state quantity obtained by filtering, and returns each estimation value to the control law and orbit determination algorithm. And an extended Kalman filter in the orbit determination algorithm calculates a state predicted value based on orbit dynamics according to the obtained estimated value, corrects the state predicted value and the navigation observed value to obtain an optimal state estimation, and returns the optimal state estimated value to a control law of parameter estimation and orbit maintenance. The control law determines the thrust to maintain the trajectory based on the estimated values obtained, and the optimal estimation of the current state. Along with the control, the estimated values of all parameters are gradually accurate, and after the estimated values are introduced into a track determination algorithm, the accuracy of the predicted values based on track dynamics is higher and higher, so that the convergence of navigation filtering errors is accelerated, and the success of track maintenance is finally ensured.
The logical relationships of the three aspects of orbit maintenance, gravitational field mapping, autonomous navigation are thus shown in fig. 5. In the parameter estimation module, the spacecraft mass, the thruster coefficient and the spherical harmonic coefficient are estimated, and all estimated values are returned to the orbit position and the extended Kalman filtering. And the extended Kalman filtering calculates a state predicted value based on orbit dynamics according to the obtained estimated value, corrects the state predicted value with the navigation observed value to obtain an optimal state estimation, and returns the optimal state estimated value to the parameter estimation and orbit maintenance module. The orbit maintenance module determines the thrust and maintains the orbit according to the obtained estimated value and the optimal estimation of the current state. Along with the control, the estimated values of various parameters are gradually accurate, and after the estimated values are introduced into the EKF, the accuracy of the predicted values based on orbit dynamics is higher and higher, so that the convergence of navigation filtering errors is accelerated, and finally the success of orbit maintenance is ensured.
Claims (5)
1. An integrated method for track maintenance, track determination and real-time mapping of a gravitational field is characterized by comprising the following steps:
s1, establishing a control law of track maintenance, a track determining algorithm and an estimation law of real-time mapping of a gravitational field, wherein the method comprises the following steps of:
s1-1, establishing an orbit dynamics model
Modeling a satellite gravitational field by adopting a spherical harmonic coefficient method, and establishing a dynamic model of the spacecraft in the satellite harmonic coefficient gravitational field, wherein modeling conditions are as follows: during real-time mapping, spin axis direction of the star is unchanged, and the spacecraft mass is ignored relative to the target celestial mass;
s1-2, establishing a control law for calculating thrust required by the spacecraft from the current position, the current speed to the expected position and the expected speed, and automatically correcting a spherical harmonic coefficient by the spacecraft to achieve an estimation law for real-time mapping of a gravitational field, wherein the method comprises the following steps of:
s1-2-1, modeling thrust errors;
s1-2-2, designing a sliding mode surface;
s1-2-3, designing an estimation law, wherein the estimation law comprises the following contents:
the adaptive estimation is realized by adopting a method of correcting parameters in real time, the perturbation acceleration caused by a nominal dynamics model relative to a real model is bounded, and the upper bound of the perturbation acceleration is set as D max ∈R 3 D is gravitational acceleration caused by model error, D is less than or equal to D max The mass m of the spacecraft and the first coefficient epsilon of the thruster are calculated 1 Second coefficient epsilon of thruster 2 And upper boundary of perturbation acceleration D max For the adaptive estimator, an estimation error is defined as:
in the above-mentioned method, the step of,error of estimated value for spacecraft quality, +.>For the spacecraft mass estimation value, m is the spacecraft mass,/->Error of upper bound estimate for perturbation acceleration, +.>For the upper bound estimate of the perturbation acceleration, D max In order to perturb the upper bound of acceleration,error of estimated value of spherical harmonic coefficient of gravitational field of star, < +.>Is the estimated value of the harmonic coefficient of the gravitational field of the star lambda j Is the gravitational field spherical harmonic coefficient of the star, lambda 1 ＝C 20 ，λ 2 ＝C 21 Sequentially numbered, the number is that,
the adaptive estimation law definition formula is:
in the above-mentioned method, the step of,the change rate of the error of the spacecraft quality estimated value is represented by q, the second coefficient of the sliding mode surface, gamma, the first estimated law parameter and m 0 ＝m 0 (t) is the upper mass limit of the spacecraft, and the condition that at any moment t, m is satisfied 0 (t)≥m(t)，/>For spacecraft quality estimation, S is the sliding mode surface, sgn is the sign function, ++>Is the difference between the true acceleration and the nominal angular velocity, +.>For estimating the rate of change of the thruster coefficient, E i For the second estimated law parameter, Ω is the third estimated law parameter, < >>For the rate of change of the perturbation acceleration upper bound estimate error, Λ D For the fourth estimated law parameter, γ, Ω, E i 、Λ D All are normal numbers and are added with->For the rate of change of the spherical harmonic coefficient estimation value error, Λ λj For the fifth estimated law parameter D λj Is the spherical harmonic coefficient lambda j The resulting perturbation acceleration is referred to as "perturbation acceleration",
obtaining each spherical harmonic coefficient lambda j The estimated value and the estimated value of the star autorotation speed omega are obtained in real time, namely the star gravitational field and the Coriolis acceleration caused by the star autorotation are obtained in real time, namely the real-time mapping of the star gravitational field is realized,
D λj about lambda j Partial derivative d of (2) λj The solution formula of (2) is:
this can be achieved by:
in the above-mentioned method, the step of,for the rate of change of the spherical harmonic coefficient estimation value error, Λ λj For the fifth estimated law parameter d λj For D λj About lambda j Q is the second coefficient of the sliding mode surface, m 0 For the upper mass boundary of the spacecraft, S is the sliding mode surface, sgn is the sign function, ++>For the spherical harmonic coefficient estimation value, < >>To estimate the thrust, u t For the theoretical thrust calculated by the control law, < +.>For the estimation of the thruster coefficients, < >>For the first thruster coefficient estimation, +.>For the second thruster coefficient estimation, +.>As the estimated value of the quality of the spacecraft,
due to the actual thrust u and the estimated thrustIs smaller and the thruster is specific to the pulse I sp Very large, it follows that:
in the above-mentioned method, the step of,g is the real mass change rate of the spacecraft 0 For the gravitational acceleration of the earth's surface, I sp For thrust specific thrust, u is the actual thrust, +.>In order to estimate the thrust force,
the rate of change of the spacecraft mass estimate can be obtained as:
in the above-mentioned method, the step of,for the rate of change of the quality estimate +.>Error rate of change for spacecraft mass estimation, < >>For the rate of change of the real mass, m 0 Is the upper mass limit of the spacecraft, S is a sliding mode surface, sgn is a sign function, g 0 Is the earth surface gravity acceleration, isp is the thruster specific impulse,>to estimate the thrust +.>Q is a second coefficient of a sliding mode surface, and gamma is a first estimation law parameter;
s1-2-4, control law design, which comprises the following contents:
the control law expression is as follows:
in the above, u t U is the theoretical thrust calculated by the control law t1 As the first term of theoretical thrust, u t2 S is a sliding mode surface, S is a theoretical thrust second term 1 Is a position error, S 2 P is a first coefficient of a sliding mode surface and is a positive constant; q is the second coefficient of the sliding die surface, and m is the normal number 0 For the upper mass limit of the spacecraft, k D For a first thrust adjustment factor, k D Not less than 1, sgn is a sign function,for the upper bound estimate of the perturbation acceleration, k ε For a second thrust adjustment factor, k ω For a third thrust adjustment factor, k s For the fourth thrust regulation factor,/- >For the first thruster coefficient estimation, +.>For the second thruster coefficient estimation, +.>Is the estimated value of the rotation speed of the star body, D λj Is the spherical harmonic coefficient lambda j Induced perturbation acceleration, ε max Is the upper limit value of the relative error epsilon of the thrust forcek ε 、k ω And k s Is positive constant for adjusting epsilon i Omega and lambda j Impact on thrust, greater k ε 、k ω And k s Will decrease epsilon i Omega and lambda j Too little effect, undermining the control law's adaptivity,
the formula for defining the sign function sgn is:
sgn ν＝diag[sgn(ν 1 ) sgn(ν 2 ) sgn(ν 3 )]
in the above formula, v is any vector, v 1 V is the first component of v 2 V is the second component of v 3 A third component of v;
s1-2-5, performing stability demonstration on the control law;
step S1-2-3, step S1-2-4 and step S1-2-5 are carried out simultaneously;
s1-3, establishing an orbit determination algorithm for continuously calculating the current position and the current speed of the spacecraft to determine the current orbit of the spacecraft;
s2, tracking the set orbit of the spacecraft by the control law, an orbit determination algorithm and an estimation law.
2. An integrated method of orbit maintenance, orbit determination and real-time mapping of the gravitational field according to claim 1, wherein said step S1-2-1 comprises the following steps:
the thrust of the thruster for the spacecraft is provided: neglecting the inertia of the system, namely considering that the control action is applied without time delay, and assuming that the thrust of the thruster has no error in the direction of the force and only has error in the magnitude, establishing a model formula of the thrust error of the thruster as follows:
In the above formula, u is the actual thrust output by the thruster, u t For theoretical thrust calculated by control law, ε is the relative error of thrust, ε 1 As the first coefficient of the thruster, epsilon 2 For the second coefficient of the thruster, m is the mass of the spacecraft, and it is assumed that in ground tests the upper bound epsilon of the relative error epsilon of the thrusts can already be known max And is also provided with
3. An integrated method of orbit maintenance, orbit determination and real-time mapping of the gravitational field according to claim 1, wherein said step S1-2-2 comprises the following steps:
according to the perturbation force caused by the nominal model attraction force and nominal model error of the spacecraft in the celestial body gravitational field and the control force of the spacecraft self-propulsion system, a kinetic equation under the nominal model is established as follows:
in the above-mentioned method, the step of,is the acceleration of the spacecraft, f is the gravitational acceleration under a nominal model, D is the gravitational acceleration caused by model errors, and +.>u is the actual thrust output by the thruster, m is the mass of the spacecraft, and the definition of the gravitational acceleration f is as follows:
f＝T -1 (ω pr t)f b
in the above, f is the gravitational acceleration under the nominal model, f b Is the projection of f in the star fixedly connected coordinate system, omega pr Is an a priori estimate of the rotational angular velocity of the star,wherein omega pr Is the rotational angular velocity of the star, t is the time,
The position control error formula of the spacecraft is defined as follows:
e 1 ＝r-χ d
defining the speed control error of the spacecraft as follows:
in the above, e 1 E is the position control error of the spacecraft 2 Is the speed control error of the spacecraft, r is the position of the spacecraft,for spacecraft speed χ d ∈R 3 For nominal position +.>At the time of the nominal speed of the vehicle,
the definition formula of the slip plane s is:
in the above, S is a slip form surface, S 1 Is a position error, S 2 For the speed error, p is the first coefficient of the sliding mode surface, the value is a positive constant, q is the second coefficient of the sliding mode surface, the value is a normal number, and the target of the control law is s=0.
4. An integrated method of orbit maintenance, orbit determination and real-time mapping of the gravitational field according to claim 1, wherein said step S1-2-5 comprises the following steps:
the function for stability demonstration was constructed as follows:
V＝V 0 +V 1
in the above formula, V is a stability proving function, V 0 For the first proof of stability function, V 1 For the second proof function of stability, m is the spacecraft mass, S is the sliding mode surface, gamma is the first estimation law parameter,k is the estimated quantity of the quality of the spacecraft ε For a second thrust adjustment factor, k ω For a third thrust adjustment factor, k s For a fourth thrust regulation factor, E i For the second evaluation law parameter,/ >For the i-th estimated value of the thruster coefficient, Ω is the third estimated law parameter, ++>Is the estimated value of the rotation angular velocity of the star, Λ λj For the fifth evaluation law parameter, ++>Is the spherical harmonic coefficient lambda j Induced gravitational acceleration perturbation, Λ D For the third evaluation law parameter, ++>For the error of the upper bound estimate of the perturbation acceleration,
as can be seen, all the terms of the stability proving function V are smaller than or equal to 0, soAnd->The filling condition of (2) is s=0, and the LaSalle invariant set theorem can know that the system converges global asymptote to s=0, namely s 1 =0 and s 2 =0, both the position error and the velocity error converge to zero.
5. An integrated method of orbit maintenance, orbit determination and real-time mapping of the gravitational field according to claim 2, wherein said step S1-3 comprises the following steps:
under the action of the gravitational field of the star and the thrust of the star, the continuous state equation of the spacecraft is as follows:
in the above formula, X (t) is a state vector and X (t) = [ X y z v ] x v y v z ] T U is the thrust of the thruster and u= [ u ] x u y u z ]M is the mass of the spacecraft, W (t) is the error of the state equation,
linearizing the continuous state equation by setting a sampling step length to obtain a discretized continuous state equation:
X k ＝F k X k-1 +B k u k +w k
in the above, X k X is the state of spacecraft at time k k-1 Is the state of the spacecraft at the moment k-1, F k For state transition matrix, B k For the control matrix at time k, u k To control vector, w k Is state noise and w k ＝[w k,x w k,y w k,z w k,vx w k,vy w k,vz ] T From errors caused by inaccurate estimated values of spherical harmonic coefficients, spacecraft mass and thruster coefficients,
the observation equation for spacecraft position is:
in the above, Z k Is the observation value of the current position of the spacecraft, v k To observe noise H k In order to observe the matrix,
the spacecraft is set to directly observe the position of the spacecraft to obtain an observation matrix H of the spacecraft k The following are provided:
assumed state noise w k And observation noise v k Is zero-mean Gaussian noise, independent of each other, and let Q k ＝cov(w k )，R k ＝cov(ν k ) The recurrence equation for the extended kalman filter is obtained as:
in the above-mentioned method, the step of,to estimate the state at time k-1 for time k, F k To act on the state transition matrix at time k-1, X k-1，k-1 To estimate the state at time k-1, u k To control vector B k To act on the control vector u k Control matrix, w k Is state noise, P k，k-1 State prediction error covariance matrix from k-1 time to k time，P k-1，k-1 As the estimated error covariance matrix at the moment k-1, Q k K is a process noise covariance matrix k For Kalman gain, H k For the observation equation, R k To observe covariance matrix of noise, X k，k As an estimated value of the state at time k, Z k P is the measurement value of the state at the moment k k，k As an estimation error covariance matrix at the k moment,
setting the prior value error to be-20%, and obtaining the spherical harmonic coefficient input value of the extended Kalman filter as follows:
in the above, lambda input For the spherical harmonic coefficient values for kalman filtering,lambda is the estimated value of the spherical harmonic coefficient p Is the prior value of the spherical harmonic coefficient.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211602711.7A CN116817932B (en) | 2022-12-13 | 2022-12-13 | Integrated method for track maintenance, track determination and real-time mapping of gravitational field |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211602711.7A CN116817932B (en) | 2022-12-13 | 2022-12-13 | Integrated method for track maintenance, track determination and real-time mapping of gravitational field |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116817932A CN116817932A (en) | 2023-09-29 |
| CN116817932B true CN116817932B (en) | 2024-03-15 |
Family
ID=88115520
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211602711.7A Active CN116817932B (en) | 2022-12-13 | 2022-12-13 | Integrated method for track maintenance, track determination and real-time mapping of gravitational field |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116817932B (en) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2948860A1 (en) * | 2015-11-20 | 2017-05-20 | Thales | Orbit transfer method for a spacecraft using a continuous or quasi-continuous thrust and embedded driving system for implementing such a method |
| CN108128484A (en) * | 2017-12-18 | 2018-06-08 | 北京理工大学 | A kind of binary-star system Maneuver strategy based on linearquadratic regulator |
| CN113998150A (en) * | 2021-11-29 | 2022-02-01 | 航天东方红卫星有限公司 | Ultra-low orbit satellite full-electric propulsion orbit maintaining system |
| CN114488806A (en) * | 2022-01-21 | 2022-05-13 | 北京航空航天大学 | Continuous thrust track maintaining method based on optimal sliding mode control |
| CN114879512A (en) * | 2022-07-06 | 2022-08-09 | 南京航空航天大学 | Spacecraft formation orbit fault-tolerant control method based on learning neural network sliding mode |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| FR3018121B1 (en) * | 2014-02-28 | 2016-02-12 | Thales Sa | METHOD FOR TRACKING A TRANSFER ORBIT OR A PHASE FOR ORKING A SPATIAL VEHICLE, IN PARTICULAR AN ELECTRICAL PROPULSION, AND APPARATUS FOR IMPLEMENTING SUCH A METHOD |
| US10180686B2 (en) * | 2016-03-17 | 2019-01-15 | Mitsubishi Electric Research Laboratories, Inc. | Concurrent station keeping, attitude control, and momentum management of spacecraft |
-
2022
- 2022-12-13 CN CN202211602711.7A patent/CN116817932B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2948860A1 (en) * | 2015-11-20 | 2017-05-20 | Thales | Orbit transfer method for a spacecraft using a continuous or quasi-continuous thrust and embedded driving system for implementing such a method |
| CN108128484A (en) * | 2017-12-18 | 2018-06-08 | 北京理工大学 | A kind of binary-star system Maneuver strategy based on linearquadratic regulator |
| CN113998150A (en) * | 2021-11-29 | 2022-02-01 | 航天东方红卫星有限公司 | Ultra-low orbit satellite full-electric propulsion orbit maintaining system |
| CN114488806A (en) * | 2022-01-21 | 2022-05-13 | 北京航空航天大学 | Continuous thrust track maintaining method based on optimal sliding mode control |
| CN114879512A (en) * | 2022-07-06 | 2022-08-09 | 南京航空航天大学 | Spacecraft formation orbit fault-tolerant control method based on learning neural network sliding mode |
Non-Patent Citations (5)
| Title |
|---|
| Feng, Zhan；Zheng, Yaru；Xu, Ming.Integrated strategy of stationkeeping, autonomous navigation, and real-time geodetical recovery of gravity fields: application into asteroid Lutetia mission.NONLINEAR DYNAMICS.2021,全文. * |
| 基于引力场不对称性的三体***轨道自主导航;王亚敏;刘银雪;蒋峻;孙煜坤;张永合;;深空探测学报(第01期);全文 * |
| 应用引力梯度测量的火星中低轨道航天器自主导航;陈培;韩锦飞;赖玉敏;孙秀聪;谭龙玉;;航天器工程(第03期);全文 * |
| 陈韬亦 ; 马鹏斌 ; 李江红 ; .基于强跟踪滤波器的机动航天器跟踪定位.无线电工程.2017,(第04期),全文. * |
| 高精度重力场下回归轨道半解析优化设计;何艳超;徐明;张润宁;李志武;;宇航学报(第05期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116817932A (en) | 2023-09-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102878995B (en) | Method for autonomously navigating geo-stationary orbit satellite | |
| CN109931955B (en) | Initial alignment method of strap-down inertial navigation system based on state-dependent lie group filtering | |
| CN108827310A (en) | A kind of star sensor secondary gyroscope online calibration method peculiar to vessel | |
| CN111189442B (en) | CEPF-based unmanned aerial vehicle multi-source navigation information state prediction method | |
| CN112325886B (en) | Spacecraft autonomous attitude determination system based on combination of gravity gradiometer and gyroscope | |
| CN110849360B (en) | Distributed relative navigation method for multi-machine collaborative formation flight | |
| Abdelrahman et al. | Sigma-point Kalman filtering for spacecraft attitude and rate estimation using magnetometer measurements | |
| Leonard et al. | Absolute orbit determination and gravity field recovery for 433 eros using satellite-to-satellite tracking | |
| de La Parra et al. | Low cost navigation system for UAV's | |
| CN116105730A (en) | Angle measurement-only optical combination navigation method based on cooperative target satellite very short arc observation | |
| Song et al. | Autonomous estimation of Allan variance coefficients of onboard fiber optic gyro | |
| Al-Jlailaty et al. | Efficient attitude estimators: A tutorial and survey | |
| CN111207773A (en) | Attitude unconstrained optimization solving method for bionic polarized light navigation | |
| CN112525204B (en) | Spacecraft inertia and solar Doppler speed combined navigation method | |
| CN111578931B (en) | High-dynamic aircraft autonomous attitude estimation method based on online rolling time domain estimation | |
| CN116817932B (en) | Integrated method for track maintenance, track determination and real-time mapping of gravitational field | |
| CN103256932A (en) | Replacement and extrapolation combined navigation method | |
| Chen et al. | High precision attitude estimation algorithm using three star trackers | |
| Liu et al. | Direction/distance/velocity measurements deeply integrated navigation for venus capture period | |
| Lou et al. | MACV/Radio integrated navigation for Mars powered descent via robust desensitized central difference Kalman filter | |
| Tan et al. | A new approach for small satellite gyroscope and star tracker fusion | |
| CN112393835A (en) | Small satellite on-orbit thrust calibration method based on extended Kalman filtering | |
| CN112393741A (en) | Air alignment method of SINS/BDS integrated navigation system based on finite time sliding mode | |
| Chen et al. | Spacecraft autonomous GPS navigation based on polytopic linear differential inclusion | |
| Luo et al. | An imu/visual odometry integrated navigation method based on measurement model optimization |
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 |