CN106482896A - A kind of contactless factor of inertia discrimination method of arbitrary shape rolling satellite - Google Patents

Abstract

The contactless factor of inertia discrimination method of satellite the invention discloses a kind of arbitrary shape rolls, choose three basic frequencies and construct the trigonometric function containing undetermined parameter as the approximate analytic solution of attitude dynamic equations, and replace system variable to characterize the state of system using these undetermined parameters, and recognized using linear minimum-variance estimation method.The attitude quaternion of rolling satellite can be measured by existing technological means, with the increase of observed quantity, estimated accuracy more and more higher.Using the value of these undetermined parameters, this method directly extrapolates the inertial parameter of rolling satellite, i.e. ratio between principal moment of inertia.Compared to the conventional constant value estimation technique, this method is not limited to the parameter estimation of axial symmetry target, and the ratio between the principal moment of inertia of the Tum bling Target of arbitrary shape all can be recognized.

Description

Non-contact inertia coefficient identification method for rolling satellite in any shape

Technical Field

The invention belongs to the field of attitude estimation and parameter identification of rolling targets, relates to a space on-orbit service technology in the field of spaceflight, and particularly relates to a non-contact inertia coefficient identification method of rolling satellites in any shapes.

Background

The objects of in-orbit servicing of spacecraft are diverse, with power-failure satellites or space debris occupying a large proportion. Due to the loss of power control, the failed satellite generally falls into a free and slow rolling state around an inertia main shaft under the action of initial angular velocity and disturbance moment, if the failed satellite needs to be subjected to close-range intersection, approach and even capture, the attitude change rule of the failed satellite needs to be pre-judged, and then a proper approach angle is selected to approach the failed satellite so as to avoid the damage of the device caused by collision.

In order to accurately predict the motion state of the rolling target for a long time in the future so as to make a more reasonable capture opportunity and a capture path, the inertial parameters of the rolling target must be accurately identified by using a non-contact method before capture. This is necessary to save the time and fuel consumed by the capturing and improve the capturing success rate.

For most failed satellites, the moment of inertia may vary due to fuel consumption or structural damage, etc. In the prior art, the attitude information of an unknown satellite can be discretely measured by using a stereoscopic vision device or a laser range finder under the condition of non-contact, but the inertial parameters of the satellite are extremely difficult to accurately identify under the condition of non-contact. The most common parameter identification method at present is to use a kalman filter to perform recursive fitting by using the moment of inertia of the kalman filter as one of state variables, but the precision is very low, and the convergence of the algorithm can be guaranteed only by giving an initial value of angular velocity measurement with sufficient precision. There are documents that propose to design a filter for identification using constant parameters, but are only applicable to a special case where the roll target is a symmetric rigid body. In practical application, the rotational inertia of the rolling satellite is measured by using a contact method after the catching, so that the task process is simplified, but additional energy consumption is increased.

Disclosure of Invention

The invention aims to provide a non-contact inertia coefficient identification method of a rolling satellite in any shape, which aims to overcome the defects in the prior art. Compared with the conventional constant value estimation method, the method is not limited to parameter estimation of the axisymmetric target, and the ratio of the main inertia moments of the rolling target in any shape can be identified.

In order to achieve the purpose, the invention adopts the following technical scheme:

a non-contact inertia coefficient identification method for rolling satellites in any shapes comprises the following steps:

1) establishing an approximate solution of an equation of motion of a rolling target quaternion lambda by an approximation method;

2) extracting a constant value parameter as a state quantity x of the system through approximate solution;

3) constructing an extended Kalman filter;

4) and obtaining an estimated value of the state quantity of the system through a Kalman filter, and estimating the inertial parameters of the rolling satellite by adopting the estimated value.

Further, the step 1) is specifically as follows:

1-1) establishing a kinematic model of the attitude quaternion of the rolling object:

wherein λ is₀,λ₁,λ₂,λ₃The body principal axis coordinate system of the rolling target is relative to the four components of the attitude quaternion, omega, of the inertia coordinate system_x,ω_y,ω_zThe angular velocity of the target is projected in a body coordinate system;

1-2) neglecting high frequency small quantity to obtain an angular velocity approximate solution of the rolling target:

wherein, ω is_xm,ω_ym,ω_zmThe maximum values of the angular velocity components in the three coordinate axis directions of the main body main axis coordinate system are respectively, t is time,is one of the constant-valued parameters to be identified,_x,_y,_zthe high-order infinitesimal quantities in the coordinate axis direction are all high-order infinitesimal quantities;

1-3) constructing an approximate solution to a quaternion differential equation describing the attitude motion using a sinusoidal function containing three characteristic frequencies, the quaternion λ ═ λ₀λ₁λ₂λ₃]^TThe specific form of the solution of (a) is:

λ＝A_4×6ξ (3)

wherein A is_4×6Is a constant matrix of 4 rows and 6 columns, and

wherein, theta is also a constant parameter to be identified;

1-4) converting a differential equation of a system into an algebraic equation, simultaneously combining an equation (1) with an equation (2) and an equation (3), merging the same frequency terms by using a trigonometric function operation rule, and obtaining an equation set consisting of 56 arithmetic equations if the equation (1) is established and coefficients before each frequency are zero according to the regularity of the trigonometric function:

wherein a is a matrix A_4×6A row vector generated after spreading into a vector by row,is composed ofTheta and omega_xm，ω_ym，ω_zmThe specific form of each item of the determined constant coefficient matrix is derived from an equation (3), the subscript represents the matrix dimension, and the equation is expanded to the following formula by considering the initial value lambda (0) of the system:

wherein the value of the constant matrix D 'is derived from D' a ═ λ (0):

the formula (4) is an overdetermined equation set with constraint number larger than variable number, a is obtained by using least square criterion,and the least square solution of theta is obtained to obtain an approximate solution of the quaternion lambda relative to the time t, namely an approximate equation of the change of the attitude quaternion along with the time, namely an approximate equation of the attitude motion of the rolling target.

Further, the step 2) is specifically as follows:

the actual observation quantity eta is a quaternion of an observation coordinate system relative to an inertia system, and the relationship between the quaternion and lambda is as follows:

η＝λομ

wherein "o" represents a quaternion multiplication operation; mu is an attitude quaternion between the observation coordinate system and the body principal axis coordinate system;

since quaternion multiplication does not involve higher-order term operation, and μ is a constant, each term η is a linear combination of terms λ, and as can be seen from equation (3), the observed quantity η also has periodic properties, namely:

η＝B_4×6ξ (5)

wherein B is a constant matrix satisfying A ξ -B ξ mu, and the matrix B is expanded by rows to obtain a row vector B_24×1And taking into account the constant parameter to be identifiedAnd θ, state quantities that together constitute a constant value of the system:

wherein the ratio of b, theta,are all constant parameters to be identified.

Further, the step 3) is specifically as follows:

the system has an observation equation of

h(x)＝η

And (3) solving a partial derivative of h relative to the state quantity x to obtain an observation sensitivity matrix:

since the state quantities of the system are all constant parameters, the state transition matrix Φ is:

Φ＝Ι_14×14(7)

wherein, I_14×14An identity matrix of 14 × 14;

the error of the state quantity caused by the disturbance moment is regarded as white noise of Gaussian distribution, and the variance matrix Q of the system process noise_kIs defined as:

Q_k＝σ²Ι_14×14(8)

wherein σ²Is the distribution variance of the process noise;

substituting the formulas (6) - (8) into a general equation of the Kalman filter to obtain a step of identifying the asymmetric rolling target parameters, and constructing the generalized Kalman filter formed by taking the parameters to be identified as state quantities.

Further, the step 4) is specifically as follows:

after an estimation value of the state quantity x is obtained through a Kalman filter, the estimation value is substituted into formula (4) and least square solution is solved, and then a matrix D and an angular velocity parameter omega are obtained_xm，ω_ymAnd omega_zmAnd further solving to obtain the inertiaEstimated value of the quantity parameter: :

wherein p is_x,p_y,p_zAnd K is a first class complete elliptic integral coefficient corresponding to the elliptic function of the angular velocity as a rotational inertia ratio parameter to be identified.

Compared with the prior art, the invention has the following beneficial technical effects:

(1) the method uses the motion law of three sine functions to approximate the attitude quaternion, is not limited in a specific solution form, and can prove that the approximation method can obtain high enough precision in all motion states, so that the method has no limitation on the shape of the target and is suitable for parameter identification of rolling targets in any asymmetric shape. (2) According to the method, constant parameters are used for replacing variable parameters to serve as state parameters of the system, so that the partial derivative of the state parameters relative to time in a nominal state is zero, and when the observed time interval is large, the error of a predicted value calculated by using numerical integration can be remarkably reduced, and the estimation accuracy of the inertia parameters is improved. (3) Since the equations have a linear form, at the characteristic frequencies θ andunder the condition that the initial value is given accurately, no requirement is made on the accuracy of other initial values, the filtering divergence phenomenon caused by too low accuracy of the initial values is avoided, and the success rate of inertia parameter estimation is improved. Numerical method using fast Fourier transformTheta andthe initial value of (2) can ensure that the precision meets the requirement.

Drawings

FIG. 1 is a flow chart for constructing an approximate solution to a quaternion kinetic equation;

FIG. 2 is a diagram of an example of observation data containing noise;

FIG. 3 is an example diagram of the use of fast Fourier transforms to transform the attitude quaternion into the frequency domain;

FIG. 4 is a diagram of an example of the identified estimated and predicted time-varying trajectories of the system motion state, where (a) is the quaternion component η corresponding to the observed coordinate system₀Is a quaternion component η corresponding to the observation coordinate system₁Is a quaternion component η corresponding to the observation coordinate system₂Is a quaternion component η corresponding to the observation coordinate system₃Is the angular velocity component ω_xIs the angular velocity component ω_yIs the angular velocity component ω_zA trajectory diagram of (a);

FIG. 5 is a diagram illustrating an example of a convergence process of relative errors of inertial parameter estimates.

Detailed Description

The invention is described in further detail below:

the invention accurately measures and calculates the inertial parameters of space debris or a failed satellite in any shape in a rolling state under the non-contact condition. The main principle is as follows: firstly, three main frequencies are selected to construct a trigonometric function containing undetermined parameters to be used as an approximate solution of the attitude motion equation of the rolling target, the set of solutions is substituted into the quaternion attitude motion equation of the rolling target, the original differential equation can be converted into an algebraic equation, and the specific value of the undetermined parameters can be obtained by solving the algebraic equation. The use of these pending parameters can characterize the state of the system in place of the system variables, i.e. the attitude quaternion of the rolling satellites is expressed as a function of these constant parameters and time. The attitude quaternion of the rolling satellite can be measured by the existing technical means, the method estimates the value of the parameter to be determined in real time by linear minimum variance estimation, and the estimation precision is higher and higher along with the increase of the observed quantity. By utilizing the values of the undetermined parameters, the method directly calculates the inertial parameters of the rolling satellite, namely the ratio between the main rotational inertia. Compared with the conventional constant value estimation method, the method is not limited to parameter estimation of the axisymmetric target, and the ratio of the main inertia moments of the rolling target in any shape can be identified.

The method specifically comprises the following steps:

the method comprises the following steps: an approximate solution of the rolling target quaternion equation of motion is established by an approximation method, as shown in fig. 1. The specific process is as follows: firstly, establishing a kinematic model of the attitude quaternion of the rolling target

Wherein λ is₀,λ₁,λ₂,λ₃Four components of attitude quaternion, omega, of principal axis coordinate system relative to inertial system as rolling target_x,ω_y,ω_zProjecting the angular speed of the target in a main shaft coordinate system of the body;

writing an approximate solution to the angular velocity of the tumbling object ignoring high frequency small quantities

Wherein ω is_xm,ω_ym,ω_zmThe maximum values of the angular velocity components in the three coordinate axis directions of the body coordinate system are respectively, t is time,is one of the constant-valued parameters to be identified,_x,_y,_zthe high-order infinitesimal quantity in the coordinate axis direction is a high-frequency term, and the long-term influence on the angular velocity is zero.

The function of quaternion changing with time can be found to have three obvious peaks in the frequency domain through frequency domain analysis, so that the approximate solution of quaternion differential equation describing attitude motion is constructed by using a sine function containing three characteristic frequencies. The solution of quaternion λ is embodied in the form of

λ＝A_4×6ξ(3)

Wherein A is_4×6Is a constant matrix of 4 rows and 6 columns, and

wherein θ isIs also a constant parameter to be identified, andthe value of (c) is determined only by the angular velocity characteristic.

And converting the differential equation of the system into an algebraic equation. The formula (2) and the formula (3) are simultaneously expressed as the formula (1), the trigonometric function operation rule is used for merging the same frequency terms, the regularity of the trigonometric function is known, if the formula (1) is established, and the coefficients before each frequency are all zero, an equation set consisting of 56 arithmetic equations is obtained

Wherein a is a matrix A_4×6A row vector generated after spreading into a vector by row,is composed ofTheta and omega_xm，ω_ym，ω_zmThe specific form of each term of the determined constant coefficient matrix is derived from equation (3), and the subscript represents the dimension of the matrix. Considering the initial value of λ (0) of the system, the above equation can be extended to

Wherein the value of the constant matrix D 'is derived from D' a ═ λ (0):

equation (4) is an overdetermined system of equations with constraint numbers greater than the variable numbers, a can be found using the least squares criterion,and the least square solution of theta is obtained to obtain an approximate solution of the quaternion lambda relative to the time t, namely an approximate equation of the change of the attitude quaternion along with the time, namely an approximate equation of the attitude motion of the rolling target. In fact, since the angular velocity of the free-tumbling object satisfies the form of Jacobi elliptic function integral, the high-frequency coefficient is small, and therefore the approximate solution determined in this way has high accuracy.

Step two: and extracting the constant value parameter as the state quantity of the system. Considering the actual observation quantity eta as a quaternion of the observation coordinate system relative to the inertial system, and the relation between the quaternion and lambda

η＝λομ

Wherein "o" represents a quaternion multiplication operation, and μ is an attitude quaternion between the observation coordinate system and the principal axis coordinate system of the body. Since quaternion multiplication does not involve higher order term operations, and μ is a constant, the η terms may all be linear combinations of the λ terms. From the form of equation (3), the observed quantity η is likewise of periodic nature, i.e. has

η＝B_4×6ξ (5)

Wherein B is a constant matrix satisfying A ξ -B ξ mu, the matrix B is expanded by rows to obtain a row vector B_24×1And taking into account the constant parameter to be identifiedAnd theta, a constant state quantity x which together form the system

b,θ，All are constant parameters to be identified;

step three: and constructing an extended Kalman filter.

The system has an observation equation of

h(x)＝η

Obtaining an observation sensitivity matrix by solving a partial derivative of h relative to the state quantity x

Since the state quantities of the system are all constant parameters, the state transition matrix phi is an identity matrix

Φ＝Ι_14×14(7)

Wherein, I_14×14An identity matrix of 14 × 14;

the error of the state quantity caused by the disturbance moment is regarded as white noise of Gaussian distribution, and the variance matrix Q of the system process noise_kIs defined as

Q_k＝σ²Ι_14×14(8)

Wherein σ²Is the distribution variance of the process noise;

substituting the formulas (6) - (8) into a general equation of a Kalman filter to obtain the step of identifying the asymmetric rolling target parameters. Similar to the case of a symmetric target, the terms of the b vector in the state quantity satisfy a linear form, without giving an accurate initial value. And theta andinvolving trigonometric functions, and therefore requiring the use of a discrete Fourier transform to determine the characteristic angular frequencies θ andthe initial value of (c).

Step four: after the estimation value of the state quantity x is obtained through a Kalman filter, the estimation value is substituted into formula (4) and least square solution is solved, and then matrix D and angular velocity parameter omega can be obtained_xm，ω_ymAnd omega_zmAn estimate of (d). Further, the estimated value of the moment of inertia parameter is calculated by the following formula

p_x,p_y,p_zAnd K is a first class complete elliptic integral coefficient corresponding to the elliptic function of the angular velocity as a rotational inertia ratio parameter to be identified.

For a better understanding of the objects and advantages of the present invention, reference is made to the following description taken in conjunction with the accompanying drawings and examples in which:

assuming that the rolling object does free-floating motion in space, the function of the four variables of the attitude quaternion along with the time can be measured by using a stereo vision device or a laser range finder, as shown in figure 2. The measurement results are contaminated by noise due to disturbing moments and observation errors. By applying the method, the inertial parameter p of the rolling satellite can be estimated in real time by using the observation noises_x，p_yAnd p_zThe method specifically comprises the following steps:

the method comprises the following steps: firstly using a fast Fourier transform algorithmλ₀Part of the data of (t) is changed to the frequency domain as shown in fig. 3. There are three peaks in the frequency domain, with the value of the abscissa of the middle peak being assigned θ, and the interval between the abscissas of the peaks being assignedAs an initial value. The initial values of the other parameters in x are assigned to 0.

Step two: and taking x as a state parameter and taking the real-time observed quantity of the attitude quaternion lambda as input, constructing a Kalman filter, and gradually estimating a more accurate value of the state parameter. As shown in fig. 4, the errors of the attitude quaternion and the angular velocity calculated and predicted from the estimated values of the constant parameters are smaller than those of the directly observed values.

Step three: calculating rolling satellite inertia parameter p by using values of x parameters_x，p_yAnd p_zThe curve of the estimated value of (a) with time is shown in fig. 5, and it can be seen that the estimated value of the inertial parameter approaches the true value as the observed quantity increases.

The values of the system parameters used in this example are shown in table 1:

TABLE 1 values of System parameters used in the examples

Claims (5)

1. A non-contact inertia coefficient identification method for a rolling satellite with any shape is characterized by comprising the following steps:

1) establishing an approximate solution of an equation of motion of a rolling target quaternion lambda by an approximation method;

2) extracting a constant value parameter as a state quantity x of the system through approximate solution;

3) constructing an extended Kalman filter;

4) and obtaining an estimated value of the state quantity of the system through a Kalman filter, and estimating the inertial parameters of the rolling satellite by adopting the estimated value.

2. The method for identifying the non-contact inertia coefficient of the rolling satellite with the arbitrary shape according to claim 1, wherein the step 1) is specifically as follows:

1-1) establishing a kinematic model of the attitude quaternion of the rolling object:

$\begin{array}{c}{\stackrel{\·}{\mathrm{\λ}}}_0=\frac{1}{2}\left(-{\mathrm{\ω}}_x{\mathrm{\λ}}_1-{\mathrm{\ω}}_y{\mathrm{\λ}}_2-{\mathrm{\ω}}_z{\mathrm{\λ}}_3\right)\\ {\stackrel{\·}{\mathrm{\λ}}}_1=\frac{1}{2}\left({\mathrm{\ω}}_x{\mathrm{\λ}}_0+{\mathrm{\ω}}_z{\mathrm{\λ}}_2-{\mathrm{\ω}}_y{\mathrm{\λ}}_3\right)\\ {\stackrel{\·}{\mathrm{\λ}}}_2=\frac{1}{2}\left({\mathrm{\ω}}_y{\mathrm{\λ}}_0-{\mathrm{\ω}}_z{\mathrm{\λ}}_1+{\mathrm{\ω}}_x{\mathrm{\λ}}_3\right)\\ {\stackrel{\·}{\mathrm{\λ}}}_3=\frac{1}{2}\left({\mathrm{\ω}}_z{\mathrm{\λ}}_0+{\mathrm{\ω}}_y{\mathrm{\λ}}_1-{\mathrm{\ω}}_x{\mathrm{\λ}}_2\right)\end{array}---\left(1\right)$

wherein λ is₀,λ₁,λ₂,λ₃The body principal axis coordinate system of the rolling target is relative to the four components of the attitude quaternion, omega, of the inertia coordinate system_x,ω_y,ω_zThe angular velocity of the target is projected in a body coordinate system;

1-2) neglecting high frequency small quantity to obtain an angular velocity approximate solution of the rolling target:

wherein, ω is_xm,ω_ym,ω_zmThe maximum values of the angular velocity components in the three coordinate axis directions of the main body main axis coordinate system are respectively, t is time,is one of the constant-valued parameters to be identified,_x,_y,_zthe high-order infinitesimal quantities in the coordinate axis direction are all high-order infinitesimal quantities;

1-3) constructing an approximate solution to a quaternion differential equation describing the attitude motion using a sinusoidal function containing three characteristic frequencies, the quaternion λ ═ λ₀λ₁λ₂λ₃]^TThe specific form of the solution of (a) is:

λ＝A_4×6ξ (3)

wherein A is_4×6Is a constant matrix of 4 rows and 6 columns, and

wherein, theta is also a constant parameter to be identified;

1-4) converting a differential equation of a system into an algebraic equation, simultaneously combining an equation (1) with an equation (2) and an equation (3), merging the same frequency terms by using a trigonometric function operation rule, and obtaining an equation set consisting of 56 arithmetic equations if the equation (1) is established and coefficients before each frequency are zero according to the regularity of the trigonometric function:

wherein a is a matrix A_4×6A row vector generated after spreading into a vector by row,is composed ofTheta and omega_xm，ω_ym，ω_zmThe specific form of each item of the determined constant coefficient matrix is derived from an equation (3), the subscript represents the matrix dimension, and the equation is expanded to the following formula by considering the initial value lambda (0) of the system:

wherein the value of the constant matrix D 'is derived from D' a ═ λ (0):

$D^{\′}=\left[\begin{array}{cccccccccccccccccccccccc}1& 0& 1& 0& 1& 0& & & & & & & & & & & & & & & & & & \\ & & & & & & 1& 0& 1& 0& 1& 0& & & & & & & & & & & & \\ & & & & & & & & & & & & 1& 0& 1& 0& 1& 0& & & & & & \\ & & & & & & & & & & & & & & & & & & 1& 0& 1& 0& 1& 0\end{array}\right]$

the formula (4) is an overdetermined equation set with constraint number larger than variable number, a is obtained by using least square criterion,and the least square solution of theta is obtained to obtain an approximate solution of the quaternion lambda relative to the time t, namely an approximate equation of the change of the attitude quaternion along with the time, namely an approximate equation of the attitude motion of the rolling target.

3. The method for identifying the non-contact inertia coefficient of the rolling satellite with the arbitrary shape according to claim 2, wherein the step 2) is specifically as follows:

the actual observation quantity eta is a quaternion of an observation coordinate system relative to an inertia system, and the relationship between the quaternion and lambda is as follows:

whereinRepresenting a quaternion multiplication; mu is an attitude quaternion between the observation coordinate system and the body principal axis coordinate system;

since quaternion multiplication does not involve higher-order term operation, and μ is a constant, each term η is a linear combination of terms λ, and as can be seen from equation (3), the observed quantity η also has periodic properties, namely:

η＝B_4×6ξ (5)

wherein B is a group ofThe matrix B is expanded according to rows to obtain row vectors B_24×1And taking into account the constant parameter to be identifiedAnd θ, state quantities that together constitute a constant value of the system:

wherein,are all constant parameters to be identified.

4. The method for identifying the non-contact inertia coefficient of the rolling satellite with the arbitrary shape according to claim 2, wherein the step 3) is specifically as follows:

the system has an observation equation of

h(x)＝η

And (3) solving a partial derivative of h relative to the state quantity x to obtain an observation sensitivity matrix:

since the state quantities of the system are all constant parameters, the state transition matrix Φ is:

Φ＝Ι_14×14(7)

wherein, I_14×14An identity matrix of 14 × 14;

the error of the state quantity caused by the disturbance moment is regarded as white noise of Gaussian distribution, and the variance matrix Q of the system process noise_kIs defined as:

Q_k＝σ²Ι_14×14(8)

wherein σ²As process noiseThe distribution variance of (c);

substituting the formulas (6) - (8) into a general equation of the Kalman filter to obtain a step of identifying the asymmetric rolling target parameters, and constructing the generalized Kalman filter formed by taking the parameters to be identified as state quantities.

5. The method for identifying the non-contact inertia coefficient of the rolling satellite with the arbitrary shape according to claim 4, wherein the step 4) is specifically as follows:

after an estimation value of the state quantity x is obtained through a Kalman filter, the estimation value is substituted into formula (4) and least square solution is solved, and then a matrix D and an angular velocity parameter omega are obtained_xm，ω_ymAnd omega_zmAnd further solving to obtain an estimated value of the inertia parameter: :

wherein p is_x,p_y,p_zAnd K is a first class complete elliptic integral coefficient corresponding to the elliptic function of the angular velocity as a rotational inertia ratio parameter to be identified.