On-Orbit Estimation Method for Spacecraft Mass Property Parameters Based on Thrusters

Abstract

The present invention relates to an on-orbit estimation method for spacecraft mass characteristic parameters based on thrusters, including: spacecraft actuator configuration, establishment of a dynamic model, establishment of a kinematic model, establishment of an attitude measurement component model, on-orbit estimation of the three-axis moment of inertia based on thrusters, on-orbit estimation of the center of mass in the steady state mode and the orbital maneuver mode, and on-orbit estimation of the spacecraft mass. The present invention solves the problems of on-orbit estimation of the spacecraft moment of inertia, center of mass, and mass, and lays a solid technical foundation for the on-orbit redesign or correction of attitude and orbit control algorithms that rely on mass characteristic parameters. Through simulation analysis, this method is reasonable and effective, with a simple algorithm, small computational load, low fuel consumption during implementation, and is easy to implement in engineering. At the same time, a flexible model is added to the simulation of this method, so it is also applicable to flexible spacecraft.

Description

Technical Field

[0001] The present invention belongs to the field of spacecraft attitude and orbit control, and particularly relates to a method for on-orbit estimation of spacecraft mass property parameters based on thrusters, which is applicable to the adaptive control of variable-structure spacecraft attitude and orbit. Background Art

[0002] With the continuous development of space technology, the structure of spacecraft has become more and more complex, undertaking more and more tasks, and the requirements for control accuracy have also become higher and higher. The complexity of the structure is reflected in that spacecraft usually carry large flexible solar panels, large flexible antennas and various payloads; the diversity of tasks leads to frequent changes in the structure of spacecraft, such as spacecraft payload release, satellite recovery, docking of space station with manned spacecraft or cargo spacecraft, etc. Generally, the design of control algorithms depends on mass property parameters, but the changes in the on-orbit structure of spacecraft make these parameters unknown. In order to be able to redesign or correct the attitude and orbit control algorithms of spacecraft to improve control accuracy and save fuel, on-orbit estimation of dynamic parameters must be carried out, and the control accuracy depends to a large extent on the accuracy of parameter estimation. Summary of the Invention

[0003] The technical problem to be solved by the present invention is: to overcome the deficiencies of the prior art and provide a method for on-orbit estimation of spacecraft mass property parameters based on thrusters.

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

[0005] A method for on-orbit estimation of spacecraft mass property parameters based on thrusters, comprising:

[0006] Taking thrusters as actuators, determining the configuration and layout of spacecraft actuators;

[0007] Establishing the rigid body dynamics equation of a spacecraft with thrusters as actuators, and considering the flexibility of the solar array, establishing the flexible dynamics equation of the spacecraft;

[0008] Establishing the kinematic equation described by Euler angles in the spacecraft body coordinate system;

[0009] Determining the configuration and layout of spacecraft attitude measurement components, and establishing the mathematical model of spacecraft attitude measurement components;

[0010] Conducting on-orbit estimation of the three-axis moment of inertia based on thrusters;

[0011] Conducting on-orbit estimation of the center of mass in the steady state mode and the orbit maneuver mode;

[0012] Conducting spacecraft mass estimation.

[0013] Preferably, in the body coordinate system, with thrusters as the actuators, the rigid body dynamics equation of the spacecraft is:

[0014]

[0015] Considering the flexibility of the solar array, the flexible dynamics equation of the spacecraft is:

[0016]

[0017] In the formula, ω is the angular velocity vector of the spacecraft body relative to the inertial space; J is the moment of inertia matrix of the spacecraft body; T c is the control torque vector; T d is the disturbance torque vector; F s is the coupling coefficient matrix of the solar array to the spacecraft; F s T is the coupling coefficient matrix of the spacecraft to the solar array; is the damping coefficient matrix of the solar array; Λ is the modal frequency matrix of the solar array; η is the modal coordinate vector of the solar array.

[0018] Preferably, the kinematic equation of the spacecraft is:

[0019] ω = ω bo + ω oi

[0020] Projected onto the spacecraft body coordinate system:

[0021]

[0022] When the attitude of the spacecraft changes by a small angle, it is simplified to:

[0023]

[0024] In the formula, ω bo represents the component of the rotational angular velocity of the spacecraft body coordinate system relative to the centroid orbit coordinate system on the spacecraft body coordinate system; ω oi represents the component of the rotational angular velocity of the centroid orbit coordinate system relative to the geocentric inertial coordinate system in the centroid orbit coordinate system; is the coordinate transformation matrix from the centroid orbit coordinate system to the body coordinate system; θ and ψ respectively represent the roll angle, pitch angle, and yaw angle of the spacecraft body coordinate system relative to the centroid orbit coordinate system, ω0 represents the orbital angular velocity, which is a constant, ω x , ω y , ω z respectively represent the components of the angular velocity vector of the spacecraft body relative to the inertial space on the three axes of the spacecraft body coordinate system.

[0025] Preferably, the implementation method for establishing the mathematical model of the spacecraft attitude measurement component is as follows:

[0026] The gyro measurement model is:

[0027]

[0028] In the formula, U is the gyro output value; ω′ is the true attitude angular velocity of the spacecraft; v g , v b is white noise with zero mean and uncorrelated with each other; b is the constant drift of the gyro;

[0029] The error δf of the accelerometer is:

[0030]

[0031] In the formula, f is the output value of the actual acceleration, Δk ai (i = x, y, z) is the error vector caused by the scale factor on the three axes x, y, z of the body coordinate system, α yx , α zx , α xy , α zy , α xz , α yz is the installation error; δb A is the random error vector. Let be the random errors of the three-axis accelerometers, then there is:

[0032]

[0033]

[0034] In the above formula, is a random constant, is a first-order Markov process, is a white noise process; α mi , are the autocorrelation parameter of the first-order Markov process and white noise.

[0035] Preferably, the on-orbit estimation method for the rolling axis moment of inertia is as follows:

[0036] Compared with the thrust moment, the environmental moment can be ignored, and the spacecraft dynamics equation becomes:

[0037]

[0038] In the formula, ρ i is the position vector of the action point of the i-th thruster relative to the centroid coordinate system, F iis the thrust vector of the i-th thruster; ω is the angular velocity vector of the spacecraft body relative to the inertial space; J is the moment of inertia matrix of the spacecraft body;

[0039] Select a suitable rolling-axis thruster to work, with no thrust acting in the pitch-axis and yaw-axis directions, so that the rolling angular velocity ω x undergoes a certain change, while the angular velocities ω y and ω z of the other two axes do not change significantly. Design the thruster jet logic. Through the measurement of the rolling angular velocity, use the least-squares method to estimate the rolling-axis moment of inertia ω x ;

[0040] The principle of on-orbit estimation of the pitch-axis and yaw-axis moments of inertia is the same as that of the rolling-axis moment of inertia on-orbit estimation.

[0041] Preferably, the on-orbit estimation method of the center of mass in the steady state mode based on thrusters is as follows:

[0042] In the steady state mode, the rigid body dynamics equation with thrusters as actuators is:

[0043]

[0044] where ω is the angular velocity vector of the spacecraft body relative to the inertial space; J is the moment of inertia matrix of the spacecraft body; let r i be the coordinates of the i-th thruster relative to the layout coordinate system, and the center of mass vector is expressed as r cm , r cm is the center of mass position vector to be estimated; F i is the thrust vector of the i-th thruster; T d is the disturbance torque vector;

[0045] Select the thrusters parallel to the O p x p axis of the spacecraft layout coordinate system to work. If J, ω, and r i and F i are known, ignoring the disturbance torque, the coordinates of the y-axis and z-axis in the center of mass layout coordinate system in the steady state mode can be estimated according to the least-squares method and

[0046] Using the same method, select the thrusters parallel to the O p y p axis to work, and the coordinates of the x-axis and z-axis in the center of mass layout coordinate system in the steady state mode can be estimated according to the least-squares method and

[0047] Preferably, the method for estimating the center of mass in the steady state mode based on "gyroscope + accelerometer" is as follows:

[0048] The output a of the accelerometer installed on the spacecraft includes two parts: translational and rotational:

[0049]

[0050] In the formula, a cm is the translational acceleration of the center of mass, r is the position vector of the thruster action point relative to the center of mass; ω is the angular velocity vector of the spacecraft body relative to the inertial space;

[0051] In the steady state mode, when there is no active force acting on the spacecraft, the translational acceleration a cm = 0; The accelerometer is installed on a rigid body, and the relative distance from the center of mass is fixed, so there is

[0052] Define the following variables:

[0053]

[0054]

[0055] ω x , ω y , ω z respectively represent the components of the angular velocity vector of the spacecraft body relative to the inertial space on the three axes of the spacecraft body coordinate system;

[0056] Then there is:

[0057] a = (H1 + H2)r

[0058] The least squares method can be applied to obtain the estimated value of the center of mass After obtaining, according to the installation position r of the thruster a The position r of the center of mass in the layout coordinate system is obtained according to the following formula cm :

[0059]

[0060] Preferably, the spacecraft mass is estimated by estimating the propellant consumption, and the method is as follows:

[0061] The propellant consumption estimation formula is as follows:

[0062]

[0063]

[0064] In the formula, F is the magnitude of the engine thrust during the orbit transfer process, Δt is the thrust action time, Ve Here, \(v\) is the exhaust velocity or specific impulse at the engine outlet, \(m_0\) is the total mass before orbit transfer, \(\Delta m\) is the fuel consumption; \(\Delta v\) is the velocity increment;

[0065] The formula for estimating the remaining propellant amount is as follows:

[0066]

[0067] In the formula, \(M\) is the remaining propellant amount in the tank, \(V\) T is the volume of the tank, \(m_0'\) is the fuel filling amount in the tank, \(P_0\) and \(T_0\) are the gas pressure and absolute temperature during filling, \(\rho_0\) is the density of the propellant during filling, \(P\) and \(T\) are the pressure and temperature of the system tank, and \(\rho\) is the density of the propellant;

[0068] Estimate the spacecraft mass based on the released payload mass, the docking spacecraft mass, and the remaining propellant amount.

[0069] Preferably, estimate the spacecraft mass through two center-of-mass estimations, and the method is as follows:

[0070] If there are large movable attachments on the spacecraft, estimate the mass through two center-of-mass estimations before and after the attachment movement;

[0071] The initial center-of-mass expression is:

[0072]

[0073] In the formula, \(M\) c is the known mass of the movable attachment, \(r\) c is the initial coordinate of the center of mass of the known movable attachment relative to the layout coordinate system, \(M\) s is the mass of the unknown other part, that is, the spacecraft mass, \(r\) s is the coordinate of the center of mass of the unknown other part relative to the layout coordinate system;

[0074] After the attachment moves, the center-of-mass expression is:

[0075]

[0076] In the formula, \(r\) c ' is the coordinate of the center of mass of the known movable attachment after movement relative to the layout coordinate system;

[0077] After measuring \(R\) and \(R'\), calculate the spacecraft mass:

[0078]

[0079] Preferably, the spacecraft layout coordinate system \(O\) p X p Y p Z p, the origin is taken at the geometric center of the spacecraft-rocket separation surface, O p Z p The axis is within the spacecraft-rocket separation surface and vertically points to the spacecraft towards the ground; O p X p The axis is perpendicular to the spacecraft-rocket separation surface and points to the payload compartment; O p Y p The axis forms a right-hand rule with the other two axes;

[0080] The body coordinate system of the spacecraft O b X b Y b Z b , the origin is at the center of mass of the spacecraft, O b X b 、O b Y b 、O b Z b The three axes are fixed to the spacecraft; O b Z b The axis vertically points to the spacecraft towards the ground; O b X b The axis is perpendicular to the spacecraft-rocket separation surface and points to the payload compartment; O b Y b Forms a right-hand rule with the other two axes; called O b X b 、O b Y b 、O b Z b The three axes are respectively the roll axis, pitch axis and yaw axis of the spacecraft.

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

[0082] By adopting the present invention, by selecting appropriate thruster working time, it can be ensured that during the estimation process, the spacecraft does not perform large-angle attitude maneuvers, has little influence on the orbit, consumes little propellant, has a simple algorithm, small calculation amount, short estimation time and high accuracy. The accuracy of on-orbit estimation of the moment of inertia, center of mass and mass can generally reach a relative error within 5%. The present invention will play an important role in solving the problem of on-orbit estimation of the mass characteristic parameters of variable-structure spacecraft. Description of the Drawings

[0083] Figure 1 is a flowchart of the on-orbit estimation method of the mass characteristic parameters of a spacecraft mainly using thrusters;

[0084] Figure 2 is a schematic diagram of the layout of the spacecraft thrusters;

[0085] Figure 3 is the algorithm flow for estimating the moment of inertia J x ;

[0086] Figure 4 Schematic diagram of the centroid position

[0087] Figure 5 Estimation algorithm flow of the centroid under the steady-state mode Specific implementation manners

[0088] The present invention will be further described in detail below by taking a certain spacecraft as an example in conjunction with the accompanying drawings.

[0089] In view of the characteristics of large changes in the moment of inertia, centroid and mass of the spacecraft, the present invention estimates the moment of inertia, centroid and mass of the spacecraft on orbit by selecting appropriate thrusters to work, laying a solid technical foundation for the on-orbit re-design or correction of the attitude and orbit control algorithm that depends on the mass characteristic parameters.

[0090] Figure 1 The following figure shows the flow chart of the method for on-orbit estimation of the mass characteristic parameters of a spacecraft mainly using thrusters, and this method is implemented through the following 7 steps:

[0091] Step 1: Determine the configuration scheme and layout scheme of the spacecraft actuator (thruster). Figure 2 Schematic diagram of the layout of the spacecraft thrusters

[0092] Step 2: Establish the rigid body dynamics model and flexible dynamics model of the spacecraft, give the dynamic equation of the spacecraft with the thruster as the actuator, and give the flexible dynamics equation of the spacecraft considering the flexibility of the solar array.

[0093] a. Rigid body dynamics model

[0094] In the body coordinate system, the rigid body dynamics equation of the spacecraft is:

[0095]

[0096] T is the total torque in the body coordinate system, and H is the angular momentum in the body coordinate system;

[0097] Using thruster control, the vector representation of the rigid body dynamics equation of the spacecraft is as follows:

[0098]

[0099] In the formula, ω is the angular velocity vector of the spacecraft body relative to the inertial space; J is the moment of inertia matrix of the spacecraft body; T c is the control torque vector; T d is the disturbance torque vector.

[0100] The disturbance torque model is as follows:

[0101] Td = T de + T dg + T dm (3)

[0102] where T de is the aerodynamic interference torque, T dg is the gravity gradient torque, and T dm is the residual magnetic torque.

[0103] b. Flexible dynamics model

[0104] Considering the flexibility of the solar array, the flexible dynamics model of the spacecraft is expressed as:

[0105]

[0106] where J is the moment of inertia matrix of the spacecraft; ω is the angular velocity vector of the spacecraft; T c is the control torque vector of the spacecraft; T d is the interference torque vector of the spacecraft; F s is the coupling coefficient matrix of the solar array to the spacecraft; F s T is the coupling coefficient matrix of the spacecraft to the solar array; is the damping coefficient matrix of the solar array; Λ is the modal frequency matrix of the solar array; η is the modal coordinate vector of the solar array.

[0107] Step 3: Based on the rigid-body dynamics equation of the spacecraft in Step 2, establish a kinematic model described by Euler angles in the spacecraft body coordinate system and simplify it.

[0108] According to the rigid-body compound motion relationship, the space rotation angular velocity vector ω of the spacecraft is equal to the rotation angular velocity vector ω b Y b Z b of the spacecraft body coordinate system OX o Y o Z o relative to the centroid orbit coordinate system OX bo Y o Z o Z o plus the transport angular velocity vector ω e X I Y I Z I of the centroid orbit coordinate system OX oi relative to the geocentric inertial coordinate system O

[0109] ω = ω bo + ω oi (5)

[0110] Project this equation onto the spacecraft body coordinate system, and we have:

[0111]

[0112] In the formula, ω oi and ω bo respectively represent the components of the angular velocity of the spacecraft body coordinate system relative to the geocentric inertial coordinate system and the centroid orbit coordinate system on the spacecraft body coordinate system; is the coordinate transformation matrix from the centroid orbit coordinate system to the body coordinate system; θ and ψ respectively represent the roll angle, pitch angle, and yaw angle of the spacecraft body coordinate system relative to the centroid orbit coordinate system, and ω0 represents the orbital angular velocity, which is a constant.

[0113] ω oi represents the component of the angular velocity of the centroid orbit coordinate system relative to the geocentric inertial coordinate system in the centroid orbit coordinate system, and is expressed as:

[0114] ω oi =[0 -ω0 0] T (7)

[0115] In the formula, ω0 is called the orbital angular velocity.

[0116] When the attitude of the spacecraft changes within a small angle (within 10 degrees), the above formula can be simplified to:

[0117]

[0118] ω x ,ω y ,ω z respectively represent the components of the angular velocity vector of the spacecraft body relative to the inertial space on the three axes of the spacecraft body coordinate system.

[0119] Step 4, determine the configuration and layout of the spacecraft attitude measurement components (such as inertial navigation combination, star sensor, infrared horizon sensor, sun sensor, magnetometer, etc.), and establish the mathematical models of the required spacecraft attitude measurement components (such as gyroscopes, accelerometers, etc.).

[0120] a. Mathematical model of gyroscope

[0121] The output of the gyroscope measurement axis can generally be expressed as:

[0122] U=ω′+Δω (9)

[0123] In the formula, U is the gyroscope output value; ω′ is the true attitude angular velocity of the spacecraft; Δω is the interference term, including gyroscope drift, measurement noise, etc., and generally the model of the interference term is taken as:

[0124] Δω=b+d+v g(10)

[0125] where b is the constant drift of the gyroscope; d is the time-correlated drift; v g is the white noise of gyroscope measurement.

[0126] The time-correlated drift d of the gyroscope is a gradually decaying process. When the spacecraft enters the steady-state control phase, the deterministic part in d is very small compared to the constant drift b. Therefore, the random part in d can be incorporated into the white noise of gyroscope measurement, and the deterministic part in d can be compensated when correcting the constant drift of the gyroscope. Thus, the gyroscope measurement model can be expressed as:

[0127]

[0128] where v g , v b is white noise with zero mean and uncorrelated with each other.

[0129] b. Mathematical model of the accelerometer table

[0130] The error of the accelerometer is:

[0131]

[0132] where f is the output value of the actual acceleration, Δk ai (i = x, y, z) is the error vector caused by the scale factor, α yx , α zx , α xy , α zy , α xz , α yz is the installation error; δb A is the random error vector. Let be the random errors of the three-axis accelerometers, then there is:

[0133]

[0134]

[0135] In the above formula, is a random constant, is a first-order Markov process, is a white noise process. α mi , are the autocorrelation parameter of the first-order Markov process and white noise.

[0136] Step 5: Based on the rigid body dynamics equation with thrusters as the actuators in Step 2, select appropriate thrusters for the roll axis, pitch axis, and yaw axis to work, so that the angular velocity of the corresponding axis changes in a certain way, while the angular velocities of the other two axes are small quantities. By measuring the angular velocity and combining with the dynamics model, estimate the inertia of the three axes respectively.

[0137] The criteria for on-orbit estimation of the moment of inertia using thrusters are: minimizing the propellant consumption as much as possible; not causing large-angle maneuvers of the attitude; minimizing the impact on the orbit as much as possible; and having a simple algorithm and small computational load.

[0138] Since the products of inertia of the spacecraft are much smaller than the principal inertia and are easily submerged by noise if estimated, and the design of the attitude control algorithm mainly depends on the principal inertia, the products of inertia do not need to be estimated.

[0139] The dynamic equation of a rigid body with thrusters as the actuators is:

[0140]

[0141] where ρ i is the position vector of the action point of the i-th thruster relative to the centroid coordinate system, F i is the thrust vector of the i-th thruster, T d is the external disturbance torque vector, ρ i × is the cross product operation of ρ i , ω is the angular velocity matrix of the body relative to the inertial space, is the principal inertia matrix of the spacecraft.

[0142] Compared with the thrust torque, the environmental torque is very small, and T d can be ignored. The dynamic equation of the rigid body becomes:

[0143]

[0144] a. On-orbit estimation of the roll axis moment of inertia

[0145] Select appropriate thrusters for the roll axis to work, so that the roll angular velocity ω x changes in a certain way, while the angular velocities ω y and ω z of the other two axes are small quantities. The rigid body dynamics equation is written as:

[0146]

[0147] By measuring the roll angular velocity, the inertia of the roll axis can be estimated. Considering that there is no thrust in the y-axis and z-axis directions and only environmental disturbance torques exist, ω y and ω z are small quantities, J xy 、Jxz is a small quantity compared with J x , and its corresponding term can be treated as white noise. The rigid body dynamics equation can be:

[0148]

[0149] For Figure 2 the spacecraft shown, thrusters 9a and 15a are selected to work, generating a rolling acceleration in the +X direction. When thrusters 9a and 15a are working, there are:

[0150]

[0151] where r i is the coordinate of the i-th thruster relative to the layout coordinate system, and r cm is the coordinate of the center of mass relative to the layout coordinate system. The thrusts of each thruster are the same, and there are:

[0152] F 15 =-F9(20)

[0153]

[0154] Define the following variables:

[0155]

[0156] θ k =J x (23)

[0157] y k =ω T [(r9 - r 15 )×F9](24)

[0158] Then there is:

[0159] y k =φ k θ k +v k (25)

[0160] y k is the measured value, θ k is the estimated value, v k is the noise, and φ k is the coefficient.

[0161] Take the performance index function J f :

[0162]

[0163] Obtain the estimated value that minimizes J f ​is:

[0164]

[0165] In the formula, Φ = [φ1, φ2, …, φ N T , Y = [y1, y2, …, y N T .

[0166] Φ is an array composed of φ k , that is, an array composed of ω x (t) , ω x (t) is measurable, is not measurable. When the thruster is working, within one attitude control cycle, the attitude angular velocity is relatively large and the star sensor is unavailable. At this time, the estimated value of the angular velocity of the gyroscope cannot be obtained from the attitude determination system, but can be obtained through the following method:

[0167] U x = ω x + b x + v gx (28)

[0168] U x is the measured value of the angular velocity of the x-axis in the body coordinate system;

[0169] The estimated value of the angular velocity of the gyroscope is:

[0170]

[0171] In the formula, can be obtained from the attitude determination system in the steady state mode.

[0172] is not measurable and is constructed by the following method:

[0173]

[0174] During the process of estimating the moment of inertia, the design of the thruster jet logic is very important. It is necessary to avoid large-angle attitude maneuvers of the spacecraft, avoid gyro saturation, and not exceed the gyro measurement range. For the Figure 2 shown spacecraft, if two thrusters are fully jetting, the thrusters can only work for one control cycle, but it is difficult to effectively estimate J x using only the data of one control cycle. Therefore, it is necessary to design the command jet time of the thruster, and the method is as follows:

[0175]

[0176] In the formula, ω max ​​is the maximum value of the given attitude angular velocity.

[0177] Y is an array composed of y k and can be calculated from the angular velocity, the position of the thruster action point, and the thrust value. The attitude angular velocity is measurable, the coordinates of the thruster action point relative to the layout coordinate system are known, and the magnitude of the thrust can be calibrated in orbit. Therefore, y k can be obtained. k

[0178] For Figure 2 the spacecraft shown, the algorithm flow for estimating the moment of inertia J x is as shown in Figure 3 and the specific steps are as follows:

[0179] First step, give the thruster start working instruction and input the initial values t0, ω x0 , b x , r9, r 15 ;

[0180] Second step, within the first and second control cycles of the simulation, calculate φ1 and φ2 according to Equation (30), calculate y1 and y2 according to Equation (22);

[0181] Third step, within the third control cycle of the simulation, calculate φ3 according to Equation (30), calculate y3 according to Equation (22);

[0182] Fourth step, within the fourth and fifth control cycles of the simulation, calculate φ4 and φ5 according to Equation (30), calculate y4 and y5 according to Equation (22);

[0183] Fifth step, obtain five groups of data, and calculate the estimated value of the moment of inertia about the x-axis according to Equation (27)

[0184] b. On-orbit estimation of the moment of inertia about the pitch axis

[0185] Adopt a method similar to the on-orbit estimation of the moment of inertia about the roll axis. For Figure 2 the spacecraft shown, select thrusters 12a and 15a to work, generate a pitch acceleration along the +Y direction, and the rigid body dynamics equation can be written as:

[0186]

[0187] Define the following variables:

[0188]

[0189] θ​k = J y (34)

[0190]

[0191] Then we have:

[0192] y k = φ k θ k + v k (36)

[0193] y k is the measured value, θ k is the estimated value, v k is the noise, φ k is the coefficient.

[0194] The measurement equation of the pitch axis of the gyroscope is:

[0195] U y = ω y + b y + v gy (37)

[0196] U y is the measured value of the angular velocity of the y-axis in the body coordinate system;

[0197] The estimated value of the angular velocity of the gyroscope is:

[0198]

[0199] In the formula, can be obtained from the attitude determination system in the steady state mode.

[0200] is unmeasurable and is constructed by the following method:

[0201]

[0202] The design of the thruster jet logic, the calculation of y k and the calculation of the estimated value are similar to the method of the roll axis.

[0203] c. On-orbit estimation of the yaw axis moment of inertia

[0204] Adopt a method similar to the on-orbit estimation of the roll axis and pitch axis moments of inertia. For the Figure 2 spacecraft shown in, select thrusters 6a and 7a to work, generate a yaw acceleration along the +Z direction, and the rigid body dynamics equation can be written as:

[0205]

[0206] Define the following variables:

[0207]

[0208] θ k = J z (42)

[0209]

[0210] Then we have:

[0211] y k = φ k θ k + v k (44)

[0212] y k is the measured value, θ k is the estimated value, v k is the noise, and φ k is the coefficient.

[0213] U z is the measured value of the angular velocity of the y-axis in the body coordinate system;

[0214] The measurement equation of the pitch axis of the gyroscope is:

[0215] U z = ω z + b z + v gz (45)

[0216] The estimated value of the angular velocity of the gyroscope is:

[0217]

[0218] In the formula, can be obtained from the attitude determination system in the steady state mode.

[0219] is unmeasurable and is constructed by the following method:

[0220]

[0221] The design of the thruster jet logic, the calculation of y k and the calculation of the estimated value are similar to the methods for the roll axis and the pitch axis.

[0222] Step 6: Based on the rigid body dynamics equation with thrusters as the actuators in Step 2, by selecting appropriate thruster operations, the on-orbit estimation of the center of mass can be carried out in the steady state mode and the orbit maneuver mode according to the least squares method. Meanwhile, the on-orbit estimation of the center of mass in the steady state mode can also be performed in the way based on "gyro + accelerometer" according to the least squares method.

[0223] a. On-orbit estimation of the center of mass based on thrusters in the steady state mode

[0224] For Figure 2 the shown spacecraft, the schematic diagram of the center of mass position is as Figure 4 shown. The dynamic equation of the rigid body with thrusters as the actuators is:

[0225]

[0226] In the formula, let r i be the coordinate of the i-th thruster relative to the layout coordinate system, and the center of mass vector is expressed as r p x p y p z p in the layout coordinate system O cm = [x c y c z c T where r cm is the center of mass position vector to be estimated. The above equation can be rewritten as:

[0227]

[0228] To obtain y c and z c , thrusters parallel to O p x p need to work. Here, thrusters 1a and 2a are selected. If J, ω and can be obtained, r i and F i are known, and the disturbance torque is ignored, then y c and z c can be estimated. The above dynamic equation can be written as:

[0229]

[0230] Define the following variables:

[0231] φ k = (F 1a + F 2a ) × (51)

[0232] θ k = r cm ​(52)

[0233]

[0234] Then there is:

[0235] y k = φ k θ k + v k (54)

[0236] y k is the measured value, θ k is the estimated value, v k is the noise, φ k is the coefficient.

[0237] According to the least squares method, and can be estimated. Since the thrusts of thrusters 1a and 2a are parallel to the O p x p axis, therefore, for the centroid offset x p x p in the direction of the O c axis, it cannot be estimated.

[0238] Using the same method, select thrusters 5a and 6a that are parallel to the O p y p to work. According to the least squares method, and

[0239] When 1a, 2a or 5a, 6a work simultaneously, the influence on the attitude is small. Therefore, the attitude angular velocity can be obtained through the attitude determination system, reducing the influence of gyro noise on the estimation accuracy.

[0240] The design of the thruster jet logic is similar to the design method in the moment of inertia estimation.

[0241] For the Figure 2 shown spacecraft, the algorithm flow for estimating the centroid in the steady state mode is as Figure 5 shown. The specific steps are as follows:

[0242] First step, give the thruster start working instruction and input the initial values t0, ω0, r 1a , r 2a ;

[0243] Second step, in the first and second control cycles of the simulation, calculate according to Equation (30) Calculate φ1 and φ2 according to Equation (51), and calculate y1 and y2 according to Equation (53);

[0244] In the third step, within the third control cycle of the simulation, calculate according to Equation (30) Calculate φ3 according to Equation (51) and calculate y3 according to Equation (53);

[0245] In the fourth step, calculate according to Equation (30) within the fourth and fifth control cycles of the simulation Calculate φ4 and φ5 according to Equation (51) and calculate y4 and y5 according to Equation (53);

[0246] In the fifth step, obtain five groups of data and calculate y c and z c estimated values of and

[0247] b. On-orbit Estimation of the Center of Mass Based on "Gyro + Accelerometer" in Steady State Mode

[0248] The output of the accelerometer installed on the spacecraft contains two parts: translational and rotational:

[0249]

[0250] where a cm is the translational acceleration of the center of mass, and r is the position vector of the thruster action point relative to the center of mass. When there is no active force acting on the spacecraft, the translational acceleration a cm = 0; the accelerometer is generally installed on the rigid body, and the relative distance from the center of mass is fixed, so there is

[0251] Define the following variables:

[0252]

[0253]

[0254] Then there is:

[0255] a = (H1 + H2)r (58)

[0256] The estimated value of the center of mass can be obtained by applying the least squares method After obtaining it, according to the installation position r of the thruster a Calculate the position r of the center of mass in the layout coordinate system according to the following formula cm :

[0257]

[0258] is unmeasurable and is still constructed by the following method:

[0259]

[0260] c. On-orbit Estimation of the Center of Mass in Orbit Maneuver Mode

[0261] The dynamic equation of a rigid body with thrusters as actuators is as follows:

[0262]

[0263] In the formula, ρ i is the position vector of the action point of the i-th thruster relative to the center-of-mass coordinate system, F i is the thrust vector of the i-th thruster, T d is the external disturbance torque vector, [ρ i ×] is the skew-symmetric matrix of ρ i , ω is the angular velocity matrix of the body relative to the inertial space, is the principal inertia matrix of the spacecraft. Let r i be the coordinates of the action point of the i-th thruster relative to the layout coordinate system, and r cm be the coordinates of the center of mass relative to the layout coordinate system, then we have:

[0264] ρ i = r i - r cm (62)

[0265] During orbit maneuver, compared with the control torque, the environmental torque is very small, and T d can be ignored, then we have:

[0266]

[0267] Define the following variables:

[0268]

[0269] θ k = r cm (65)

[0270]

[0271] y k is the measured value, θ k is the estimated value, v k is the noise, and φ k is the coefficient.

[0272] The least squares method can be used to estimate

[0273] Step 7: Perform spacecraft mass estimation. The spacecraft mass estimation can be carried out by estimating the propellant consumption or by estimating the center-of-mass positions at two times.

[0274] a. Estimating the center of mass through the estimation of propellant consumption

[0275] The reasons affecting the mass change of the spacecraft mainly include payload release, docking of the space station complex, and fuel consumption of the propellant. Only the propellant consumption is unknown. By estimating the propellant consumption, the on-orbit estimation of the spacecraft mass can be completed.

[0276] The method for estimating the propellant consumption is as follows:

[0277]

[0278]

[0279] In the formula, F is the magnitude of the engine thrust during the orbit transfer process, Δt is the thrust action time, V e is the exhaust velocity or specific impulse at the engine outlet, m0 is the total mass before the orbit transfer, Δm is the fuel consumption, and Δv is the velocity increment.

[0280] The method for estimating the remaining propellant is as follows:

[0281]

[0282] In the formula, M is the remaining amount of the propellant in the storage tank, V T is the volume of the storage tank, m0′ is the fuel filling amount of the storage tank, P0 and T0 are the gas pressure and absolute temperature during filling, ρ0 is the density of the propellant during filling, P and T are the pressure and temperature of the system storage tank, and ρ is the density of the propellant.

[0283] b. Estimating the mass through two center-of-mass estimations

[0284] If the spacecraft is equipped with large movable attachments, the mass can be estimated by performing two center-of-mass estimations before and after the attachment movement.

[0285] The initial center-of-mass expression is:

[0286]

[0287] In the formula, M c is the known mass of the movable attachment, r c is the initial coordinate of the center of mass of the known movable attachment relative to the layout coordinate system, M s is the unknown mass of the other parts, r s is the coordinate of the center of mass of the unknown other parts relative to the layout coordinate system.

[0288] After the attachment moves, the center-of-mass expression is:

[0289]

[0290] where r' c is the coordinate of the centroid of the known movable attachment after movement relative to the layout coordinate system.

[0291] After measuring R and R', the mass of the spacecraft can be calculated:

[0292]

[0293] So far, the in-orbit estimation of the spacecraft mass characteristics is completed.

[0294] The following takes Figure 2 the spacecraft shown as an example for verification through numerical simulation:

[0295] The initial conditions of the simulation are as follows:

[0296] Moment of inertia about the roll axis: 973.2 kg·m 2 ;

[0297] Moment of inertia about the pitch axis: 969.8 kg·m 2 ;

[0298] Moment of inertia about the yaw axis: 572.4 kg·m 2 ;

[0299] The coordinates of the thruster relative to the layout coordinate system are shown in Table 1.

[0300] Table 1 Thruster installation position

[0301]

[0302] Thruster thrust magnitude: 25 N, thrust error 10%;

[0303] Gyro constant drift (3σ): 5° / h;

[0304] Gyro random drift (3σ): 0.15° / h;

[0305] Coordinates of the centroid relative to the layout coordinate system: x c = 576.44, y c = 18.36, z c = -54.88;

[0306] Coordinates of the accelerometer relative to the layout coordinate system: x a = 176.44, y a = 318.36, z a = 5;

[0307] Gyro noise characteristics: zero position 3.5° / h, scale factor 0.01° / s, installation skew 30″, white noise mean square deviation 0.004° / h, white noise mean square deviation (driven random walk) 0.15° / h;

[0308] Additive table noise characteristics: zero position 5×10 -8 g0, scale factor 1×10 -8 g0, installation skew 30″, white noise mean square deviation 5×10 -8 g0, white noise mean square deviation (driving first-order Markov process) 1×10 -8 g0;

[0309] In the orbit maneuver mode, the initial three-axis attitude: θ = 0.01°, ψ = 0.01°;

[0310] In the orbit maneuver mode, the initial three-axis attitude angular velocity:

[0311] In the orbit maneuver mode, Δa = 100 km;

[0312] In the orbit maneuver mode, Δu = 2°;

[0313] In the orbit maneuver mode, Δi = 0.41°.

[0314] The simulation results are as follows:

[0315] Under the same initial conditions, after 5 simulations, the estimated results of the rolling-axis moment of inertia are shown in Table 2, the estimated results of the pitching-axis moment of inertia are shown in Table 3, and the estimated results of the yaw-axis moment of inertia are shown in Table 4.

[0316] Table 2 Estimated results of the rolling-axis moment of inertia

[0317]

[0318] Table 3 Estimated results of the pitching-axis moment of inertia

[0319]

[0320] Table 4 Estimated results of the yaw-axis moment of inertia

[0321]

[0322] In the steady-state mode, based on the centroid estimation of the thrusters, under the same initial conditions, after 5 simulations, the estimated results of the centroid on the Y-axis are shown in Table 5, the estimated results of the centroid on the Z-axis are shown in Table 6, and the estimated results of the centroid on the X-axis are shown in Table 7.

[0323] Table 5 Estimated results of the centroid on the Y-axis based on the thrusters in the steady-state mode

[0324]

[0325] Table 6 Thrust-based Z-axis centroid estimation results in steady state mode

[0326]

[0327] Table 7 Thrust-based X-axis centroid estimation results in steady state mode

[0328]

[0329] In steady state mode, for the centroid estimation based on "gyro + accelerometer", the estimation result is the difference between the centroid of the plus table and the centroid of the spacecraft. After 5 simulations under the same initial conditions, the X-axis centroid difference estimation results are shown in Table 8, the Y-axis centroid difference estimation results are shown in Table 9, and the Z-axis centroid difference estimation results are shown in Table 10.

[0330] Table 8 Sensitive device-based X-axis centroid difference estimation results in steady state mode

[0331]

[0332] Table 9 Sensitive device-based Y-axis centroid difference estimation results in steady state mode

[0333]

[0334] Table 10 Sensitive device-based Z-axis centroid difference estimation results in steady state mode

[0335]

[0336] In orbit maneuver mode, for the centroid estimation based on the thrust, after 5 simulations under the same initial conditions, the Y-axis centroid estimation results are shown in Table 11, the Z-axis centroid estimation results are shown in Table 12, and the X-axis centroid estimation results are shown in Table 13.

[0337] Table 11 Thrust-based Y-axis centroid estimation results in orbit maneuver mode

[0338]

[0339] Table 12 Thrust-based Z-axis centroid estimation results in orbit maneuver mode

[0340]

[0341] Table 13 Thrust-based X-axis centroid estimation results in orbit maneuver mode

[0342]

[0343] The centroid estimation algorithm based on the thrust has a relatively high estimation accuracy, while the centroid estimation algorithm based on "gyro + accelerometer" requires the zero bias and measurement noise of the accelerometer to be within 10-8 In the order of magnitude of g0, the advantages of the two methods are compared as shown in Table 14.

[0344] Table 14 Comparison of centroid estimation algorithms

[0345]

[0346] As can be seen from the above, by adopting the present invention and selecting appropriate thruster working time, it can be ensured that during the estimation process, the spacecraft does not perform large-angle attitude maneuvers, has little impact on the orbit, consumes less propellant, has a simple algorithm, small computational load, short estimation time, and high accuracy. Generally, the accuracy of on-orbit estimation of the moment of inertia, centroid, and mass can reach a relative error within 5%.

[0347] The present invention solves the problem of on-orbit estimation of the moment of inertia, centroid, and mass of a spacecraft, laying a solid technical foundation for the on-orbit redesign or correction of attitude-orbit control algorithms relying on mass characteristic parameters. Through simulation analysis, the method is reasonable and effective, has a simple algorithm, small computational load, consumes less fuel during implementation, is easy to implement in engineering. At the same time, the flexible model is added to the simulation of this method, so it is also applicable to flexible spacecraft.

Claims

1. An on-orbit estimation method for the mass characteristic parameters of a spacecraft based on thrusters, characterized in that, Including: Taking thrusters as the actuators, determine the configuration and layout of the spacecraft actuators; Establish the rigid-body dynamics equations of the spacecraft with thrusters as the actuators, and considering the flexibility of the solar array, establish the flexible dynamics equations of the spacecraft; Establish the kinematic equations described by Euler angles in the spacecraft body coordinate system; Determine the configuration and layout of the spacecraft attitude measurement components, and establish the mathematical model of the spacecraft attitude measurement components; Conduct on-orbit estimation of the three-axis moment of inertia based on thrusters; Conduct on-orbit estimation of the center of mass in the layout coordinate system in the steady-state mode and orbit maneuver mode; Conduct spacecraft mass estimation; The method for on-orbit estimation of the center of mass in the steady-state mode based on thrusters is as follows: In the steady-state mode, the rigid-body dynamics equations with thrusters as the actuators are: where ω is the angular velocity vector of the spacecraft body relative to the inertial space; J is the inertia matrix of the spacecraft body; let r i be the coordinates of the i-th thruster relative to the layout coordinate system, and the centroid vector is expressed as r cm , r cm is the centroid position vector to be estimated; F i is the thrust vector of the i-th thruster; T d is the disturbance torque vector; Select thrusters that work parallel to the spacecraft layout coordinate system O p x p If J, ω, and r i and F i are known and disturbance torques are ignored, the coordinates of the y-axis and z-axis in the centroid layout coordinate system under the steady-state mode can be estimated according to the least squares method and Using the same method, select the thrusters parallel to O p y p to work. According to the least squares method, the coordinates of the x-axis and z-axis in the centroid layout coordinate system under the steady-state mode can be estimated and 2. The on-orbit estimation method for spacecraft mass characteristic parameters based on thrusters according to claim 1, characterized in that In the body coordinate system, with thrusters as the actuators, the rigid-body dynamics equations of the spacecraft are: Considering the flexibility of the solar array, the flexible dynamics equations of the spacecraft are: where ω is the angular velocity vector of the spacecraft body relative to the inertial space; J is the moment of inertia matrix of the spacecraft body; T c is the control torque vector; T d is the disturbance torque vector; F s is the coupling coefficient matrix of the solar array to the spacecraft; F s T is the coupling coefficient matrix of the spacecraft to the solar array; is the damping coefficient matrix of the solar array; Λ is the modal frequency matrix of the solar array; η is the modal coordinate vector of the solar array.

3. The on-orbit estimation method for spacecraft mass characteristic parameters based on thrusters according to claim 1, characterized in that, The spacecraft kinematic equations are: ω = ω bo + ω oi Projected onto the spacecraft body coordinate system: When the attitude of the spacecraft changes by a small angle, it is simplified to: where ω bo represents the components of the angular velocity of the spacecraft body coordinate system relative to the centroid orbit coordinate system on the spacecraft body coordinate system; ω oi represents the components of the angular velocity of the centroid orbit coordinate system relative to the geocentric inertial coordinate system in the centroid orbit coordinate system; is the coordinate transformation matrix from the centroid orbit coordinate system to the body coordinate system; θ and ψ respectively represent the roll angle, pitch angle, and yaw angle of the spacecraft body coordinate system relative to the centroid orbit coordinate system, ω0 represents the orbital angular velocity, which is a constant, ω x , ω y , ω z respectively represent the components of the angular velocity vector of the spacecraft body relative to the inertial space on the three axes of the spacecraft body coordinate system.

4. A method for on-orbit estimation of spacecraft mass property parameters based on thrusters according to claim 1, characterized in that The implementation method for establishing the mathematical model of the spacecraft attitude measurement components is as follows: The gyro measurement model is: Wherein, U is the gyro output value; ω′ is the true attitude angular velocity of the spacecraft; v g , v b is white noise with zero mean, and they are uncorrelated with each other; b is the constant drift of the gyro The error δf of the accelerometer is: where f is the output value of the actual acceleration, and Δk ai (i = x, y, z) is the error vector caused by the scale factor in the three axes of the body coordinate system x, y, z, and α yx , α zx , α xy , α zy , α xz , α yz is the installation error; δb A is the random error vector. Let be the random errors of the three-axis accelerometers. Then we have: In the above formula, is a random constant, is a first-order Markov process, is a white noise process; α mi , are the autocorrelation parameter of the first-order Markov process and white noise.

5. The on-orbit estimation method for spacecraft mass characteristic parameters based on thrusters according to claim 1, characterized in that The method for on-orbit estimation of the moment of inertia about the roll axis is as follows: Compared with the thrust moment, the environmental moment can be ignored, and the spacecraft dynamics equations become: where ρ i is the position vector of the action point of the i-th thruster relative to the centroid coordinate system, and F i is the thrust vector of the i-th thruster; ω is the angular velocity vector of the spacecraft body relative to the inertial space; J is the inertia matrix of the spacecraft body; Select a suitable rolling axis thruster to operate. There is no thrust in the pitch axis and yaw axis directions, causing the rolling angular velocity ω x to change in some way, while the angular velocities ω y and ω z of the other two axes do not change significantly. Design the thruster jet logic. Through the measurement of the rolling angular velocity, use the least squares method to estimate the rolling axis moment of inertia ω x ; The principle for on-orbit estimation of the moment of inertia about the pitch axis and yaw axis is the same as that for on-orbit estimation of the moment of inertia about the roll axis.

6. The on-orbit estimation method for spacecraft mass characteristic parameters based on thrusters according to claim 1, characterized in that The method for on-orbit estimation of the center of mass in the steady-state mode based on "gyro + accelerometer" is as follows: The output a of the accelerometers installed on the spacecraft includes two parts: translational and rotational: where a cm is the translational acceleration of the center of mass, r is the position vector of the thruster application point relative to the center of mass; ω is the angular velocity vector of the spacecraft body relative to the inertial space; In the steady state mode, when there is no active force acting on the spacecraft, the translational acceleration a cm = 0; The accelerometer is installed on the rigid body, and the relative distance from the center of mass is fixed. There is Define the following variables: ω x , ω y , ω z respectively represent the components of the angular velocity vector of the spacecraft body relative to the inertial space on the three axes of the spacecraft body coordinate system; Then there is: a = (H1 + H2)r The estimated value of the centroid can be obtained by using the least squares method. After obtaining it, according to the installation position r of the thruster a The position r of the centroid in the layout coordinate system is calculated by the following formula cm :

7. A method for on-orbit estimation of spacecraft mass property parameters based on thrusters according to claim 1, characterized in that, The spacecraft mass is estimated by estimating the propellant consumption, and the method is as follows: The formula for estimating the propellant consumption is as follows: Where F is the magnitude of the engine thrust during the orbit transfer process, Δt is the thrust application time, V e is the exhaust velocity or specific impulse at the engine outlet, m0 is the total mass before the orbit transfer, and Δm is the fuel consumption; Δv is the velocity increment; The formula for estimating the remaining propellant is as follows: Where M is the remaining amount of propellant in the storage tank, V T is the volume of the storage tank, m0' is the fuel filling amount of the storage tank, P0 and T0 are the gas pressure and absolute temperature during filling, ρ0 is the density of the propellant during filling, P and T are the pressure and temperature of the system storage tank, and ρ is the density of the propellant; Estimate the spacecraft mass based on the released payload mass, the mass of the docking vehicle, and the remaining propellant.

8. A method for on-orbit estimation of spacecraft mass property parameters based on thrusters according to claim 1, characterized in that The spacecraft mass is estimated by conducting two center-of-mass estimations, and the method is as follows: If there are large movable attachments on the spacecraft, estimate the mass by conducting two center-of-mass estimations before and after the attachment moves; The initial center-of-mass expression is: where, M c is the known mass of the movable attachment, r c is the initial coordinate of the centroid of the known movable attachment relative to the layout coordinate system, M s is the mass of the unknown other part, i.e., the spacecraft mass, r s is the coordinate of the centroid of the unknown other part relative to the layout coordinate system; After the attachment moves, the center-of-mass expression is: where r′ c is the coordinate of the centroid of the known movable attachment after movement relative to the layout coordinate system; After measuring R and R', calculate the spacecraft mass:

9. The on-orbit estimation method for spacecraft mass characteristic parameters based on thrusters according to claim 1, characterized in that: Spacecraft layout coordinate system O p X p Y p Z p , the origin is taken at the geometric center of the spacecraft-rocket separation surface, O p Z p axis is in the spacecraft-rocket separation surface and perpendicular to the spacecraft pointing to the ground; O p X p axis is perpendicular to the spacecraft-rocket separation surface and points to the payload bay; O p Y p axis follows the right-hand rule with the other two axes; Spacecraft body coordinate system O b X b Y b Z b , with the origin at the center of mass of the spacecraft, O b X b 、O b Y b 、O b Z b The three axes are fixedly connected to the spacecraft; O b Z b axis is vertically pointed at the spacecraft to the ground; O b X b axis is perpendicular to the spacecraft-rocket separation plane and points to the payload compartment; O b Y b obeys the right-hand rule with the other two axes; called O b X b 、O b Y b 、O b Z b The three axes of X, Y, and Z are respectively the roll axis, pitch axis, and yaw axis of the spacecraft.

Images