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

Abstract

The invention discloses a non-contact inertia coefficient identification method for tumbling satellites with arbitrary shapes. Three main frequencies are selected to construct trigonometric functions containing undetermined parameters as the approximate analytical solution of the attitude dynamic equation, and these undetermined parameters are used to replace the system variable representation The state of the system is identified using the linear minimum variance estimation method. The attitude quaternion of the tumbling satellite can be measured by existing technical means, and with the increase of observations, the estimation accuracy is getting higher and higher. Using the values of these undetermined parameters, this method directly calculates the inertial parameters of the rolling satellite, that is, the ratio between the main moments of inertia. Compared with the previous constant value estimation method, this method is not limited to the parameter estimation of axisymmetric targets, and can identify the ratio between the principal moments of inertia of rolling targets with arbitrary shapes.

Description

A non-contact inertial coefficient identification method for arbitrarily shaped tumbling satellites

technical field

The invention belongs to the field of attitude estimation and parameter identification of tumbling targets, relates to space on-orbit service technology in the aerospace field, and in particular to a non-contact inertia coefficient identification method for tumbling satellites with arbitrary shapes.

Background technique

Spacecraft serve a variety of objects in orbit, among which power-failed satellites or space debris occupy a large proportion. Due to the loss of power control, the failed satellite will generally fall into a state of free and slow rolling around the main axis of inertia under the action of the initial angular velocity and disturbance torque. If it is to rendezvous, approach or even capture it at close range, it must predict the law of its attitude change. Judgment, and then choose a suitable approach angle to avoid collision and damage to the device.

In order to accurately predict the motion state of the tumbling target for a long time in the future, and to formulate a more reasonable capture timing and capture path, it is necessary to use a non-contact method to accurately identify the inertial parameters of the tumbling target before capture. This is very necessary for saving the time and fuel spent in capturing and improving the success rate of capturing.

For most invalid satellites, the moment of inertia may change due to fuel consumption or structural damage. Under the existing technology, stereo vision equipment or laser range finder can be used to discretely measure the attitude information of unknown satellites without contact, but it is extremely difficult to accurately identify the inertial parameters of satellites without contact. At present, the most commonly used parameter identification method is to use the Kalman filter to recursively fit its moment of inertia as one of the state variables, but the accuracy is very low, and it is necessary to give a sufficiently accurate initial value of angular velocity measurement to ensure the convergence of the algorithm. Some literatures propose to use constant parameters to design filters for identification, but it is only applicable to the special case that the tumbling target is a symmetrical rigid body. In practical applications, the contact method is often used to measure the moment of inertia of the tumbling satellite after capture. Although this simplifies the task process, it will increase additional energy consumption.

Contents of the invention

The purpose of the present invention is to provide a non-contact inertia coefficient identification method for tumbling satellites with arbitrary shapes, so as to overcome the above-mentioned defects in the prior art. The present invention directly calculates the inertia parameters of tumbling satellites, that is, the ratio between the main moments of inertia . Compared with the previous constant value estimation method, this method is not limited to the parameter estimation of axisymmetric targets, and can identify the ratio between the principal moments of inertia of rolling targets with arbitrary shapes.

To achieve the above object, the present invention adopts the following technical solutions:

A non-contact inertia coefficient identification method for tumbling satellites with arbitrary shapes, comprising the following steps:

1) establish the approximate solution of the motion equation of the rolling target quaternion λ by an approximate method;

2) Extract the constant parameter as the state quantity x of the system through the approximate solution;

3) Construct an extended Kalman filter;

4) Obtain the estimated value of the state quantity of the system through the Kalman filter, and use the estimated value to estimate the inertial parameters of the rolling satellite.

Further, step 1) is specifically:

1-1) Establish the kinematics model of the attitude quaternion of the rolling target:

Among them, λ ₀ , λ ₁ , λ ₂ , λ ₃ are the four components of the attitude quaternion of the main axis coordinate system of the rolling target relative to the inertial coordinate system, and ω _x , ω _y , ω _z are the angular velocity of the target in the body coordinates projections in the system;

1-2) Neglecting the small amount of high frequency, the approximate solution of the angular velocity of the rolling target is obtained:

Among them, ω _xm , ω _ym , ω _zm are the maximum values of the angular velocity components in the directions of the three coordinate axes of the principal axis coordinate system of the body respectively, t is the time, is one of the constant parameters to be identified, and ε _x , ε _y , ε _z are all high-order infinitesimal quantities in the direction of the coordinate axis;

1-3) Use the sine function containing three characteristic frequencies to construct the approximate solution of the quaternion differential equation describing the attitude motion, the specific form of the solution of the quaternion λ=[λ ₀ λ ₁ λ ₂ λ ₃ ] ^T is:

λ＝A _4×6 ξ (3)

where A _4×6 is a constant matrix with 4 rows and 6 columns, and

Among them, θ is also a constant parameter to be identified;

1-4) Transform the differential equation of the system into an algebraic equation, substitute formula (2) and formula (3) into formula (1) at the same time, and use trigonometric function operation rules to combine the same frequency items. From the regularity of trigonometric functions, if Equation (1) is established, and the coefficients before each frequency must be zero, then an equation system consisting of 56 arithmetic equations is obtained:

Where a is the row vector generated after the matrix A _4×6 is expanded into a vector by row, for the reason The constant coefficient matrix determined by θ and ω _xm , ω _ym , ω _zm , the specific form of each item is derived from formula (3), and the subscript indicates the dimension of the matrix. Considering the initial value λ(0) of the system, the upper expands to:

Among them, the value of the constant matrix D' is deduced by D'a=λ(0):

Equation (4) is an overdetermined equation system with the number of constraints greater than the number of variables, using the least squares criterion to obtain a, And the least squares solution of θ, the approximate solution of the quaternion λ about time t is obtained, that is, the approximate equation of the attitude quaternion changing with time, that is, the approximate equation of the attitude motion of the rolling target.

Further, step 2) is specifically:

The actual observation value η is the quaternion of the observation coordinate system relative to the inertial system, and the relationship between it and λ is:

η=λομ

Among them, "ο" represents the quaternion multiplication operation; μ is the attitude quaternion between the observation coordinate system and the main axis coordinate system of the body;

Since the quaternion multiplication does not involve high-order term operations, and μ is a constant, the items of η are linear combinations of the items of λ. It can be seen from formula (3) that the observed quantity η also has periodic properties, that is:

η=B _4×6 ξ (5)

Where B is a constant matrix satisfying Aξ=Bξομ, expand the matrix B by row to obtain a row vector b _24×1 , and consider the constant parameter to be identified and θ, together constitute the constant state quantity of the system:

Among them, b, θ, are constant parameters to be identified.

Further, step 3) is specifically:

The observation equation of the system is

h(x)=η

Calculate the partial derivative of h relative to the state quantity x to obtain the observation sensitivity matrix:

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

Φ＝Ι _14×14 (7)

Wherein, I _14×14 is an identity matrix of 14×14;

The error of the state quantity caused by the disturbance torque is considered to be white noise of Gaussian distribution, then the variance matrix Q _k of the system process noise is defined as:

Q _k = σ ² Ι _14×14 (8)

Among them, ^σ2 is the distribution variance of the process noise;

Substituting equations (6)-(8) into the general equation of the Kalman filter is the step of identifying the parameters of the asymmetric tumbling target, and constructing a generalized Kalman filter composed of the parameters to be identified as state variables.

Further, step 4) is specifically:

After the estimated value of the state quantity x is obtained through the Kalman filter, the estimated value is substituted into formula (4) and the least square solution is obtained, that is, the estimated value of the matrix D and the angular velocity parameters ω _xm , ω _ym and ω _zm is obtained, and further Solving for an estimate of the inertia parameter:

Among them, p _x , p _y , p _z are the moment of inertia ratio parameters to be identified, and K is the first kind of complete elliptic integral coefficient corresponding to the elliptic function of angular velocity.

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

(1) This method uses three sine functions to approximate the motion law of the attitude quaternion, and is not limited to a specific solution form. It can be proved that the approximation method can achieve high enough accuracy in all motion states. There is no restriction on the shape of the target, and it is suitable for parameter identification of tumbling targets with arbitrary asymmetric shapes. (2) This method uses constant parameters instead of variable parameters as the state parameters of the system, so that the partial derivatives of the state parameters with respect to time in the nominal state are zero. When the time interval of observations is large, it can be significantly reduced. The error of the predicted value can improve the estimation accuracy of inertia parameters. (3) Since the equation has a linear form, at the characteristic frequencies θ and When the initial value of is relatively accurate, there is no requirement for the accuracy of other initial values, which avoids the phenomenon of filter divergence caused by the low accuracy of the initial value, and improves the success rate of inertia parameter estimation. Numerical methods using fast Fourier transforms to obtain θ and The initial value of , can ensure that its accuracy meets the requirements.

Description of drawings

Fig. 1 is the flow chart of constructing the approximate solution of quaternion kinetic equation;

Figure 2 is an example diagram of observation data containing noise;

Fig. 3 is an example diagram of changing attitude quaternion to frequency domain using fast Fourier transform;

Fig. 4 is an example diagram of the trajectory of the estimated value and predicted value of the identified system motion state over time, where (a) is the trajectory diagram of the quaternion component _η0 corresponding to the observation coordinate system, and (b) is the observation coordinate The locus diagram of the quaternion component η ₁ corresponding to the system, (c) is the locus diagram of the quaternion component η ₂ corresponding to the observation coordinate system, and (d) is the locus diagram of the quaternion component η ₃ corresponding to the observation coordinate system , (e) is the track figure of angular velocity component ω _x , (f) is the track figure of angular velocity component ω _y , (g) is the track figure of angular velocity component ω _z ;

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

detailed description

The present invention is described in further detail below:

The invention accurately measures and calculates the inertial parameters of the space debris of any shape in the tumbling state or the invalid satellite under the condition of non-contact. The main principle is: firstly, three main frequencies are selected to construct a trigonometric function with undetermined parameters as the approximate solution of the rolling target attitude motion equation, and this group of solutions is substituted into the rolling target's quaternion attitude dynamic equation, and the original The differential equations are transformed into algebraic equations, and the specific values of the undetermined parameters can be obtained by solving this set of algebraic equations. These undetermined parameters can be used instead of system variables to characterize the state of the system, that is, the attitude quaternion of the rolling satellite is expressed as a function of these constant parameters and time. The attitude quaternion of the tumbling satellite can be measured by existing technical means. This method estimates the value of the undetermined parameters in real time through linear minimum variance estimation, and the estimation accuracy becomes higher and higher with the increase of observations. Using the values of these undetermined parameters, this method directly calculates the inertial parameters of the rolling satellite, that is, the ratio between the main moments of inertia. Compared with the previous constant value estimation method, this method is not limited to the parameter estimation of axisymmetric targets, and can identify the ratio between the principal moments of inertia of rolling targets with arbitrary shapes.

Method of the present invention specifically comprises the following steps:

Step 1: Establish an approximate solution to the quaternion motion equation of the tumbling target by an approximate method, as shown in Figure (1). The specific process is as follows: first, establish the kinematics model of the attitude quaternion of the rolling target

Among them, λ ₀ , λ ₁ , λ ₂ , λ ₃ are the four components of the attitude quaternion of the main axis coordinate system of the rolling target relative to the inertial system, and ω _x , ω _y , ω _z are the angular velocity of the target in the main axis coordinate system of the body projection in

Ignoring the high-frequency small amount, write the approximate solution of the angular velocity of the rolling target

Where ω _xm , ω _ym , ω _zm are the maximum values of angular velocity components in the direction of the three coordinate axes of the body coordinate system, t is time, is one of the constant parameters to be identified, ε _x , ε _y , ε _z are high-order infinitesimal quantities in the direction of the coordinate axis, ε is a high-frequency term, and its long-term influence on angular velocity is zero.

Through the frequency domain analysis, it can be found that the function of quaternion changing with time has three obvious peaks in the frequency domain, so the approximate solution of the quaternion differential equation describing attitude motion is constructed by using the sinusoidal function with three characteristic frequencies. The specific form of the solution of the quaternion λ is

λ＝A _4×6 ξ(3)

where A _4×6 is a constant matrix with 4 rows and 6 columns, and

where θ and is also a constant value parameter to be identified, and The value of is determined only by the angular velocity characteristic.

Transform the system of differential equations into algebraic equations. Substitute Equation (2) and Equation (3) into Equation (1) at the same time, and use trigonometric function operation rules to combine the same frequency items. From the regularity of trigonometric functions, it can be seen that if Equation (1) is true, the coefficients before each frequency must be are all zero, an equation system consisting of 56 arithmetic equations is obtained

Where a is the row vector generated after the matrix A _4×6 is expanded into a vector by row, for the reason The constant coefficient matrix determined by θ and ω _xm , ω _ym , ω _zm , the specific form of its items can be derived from formula (3), and the subscript indicates the dimension of the matrix. Considering the initial value λ(0) of the system, the above formula can be extended to

Among them, the value of the constant matrix D' is deduced by D'a=λ(0):

Equation (4) is an overdetermined equation system in which the number of constraints is greater than the number of variables, and a can be obtained by using the least squares criterion, And the least squares solution of θ, the approximate solution of the quaternion λ about time t is obtained, that is, the approximate equation of the attitude quaternion changing with time, that is, the approximate equation of the attitude motion of the rolling target. In fact, since the angular velocity of the free-rolling target satisfies the integral form of the Jacobi elliptic function, the coefficient of the high-frequency term is very small, so the approximate solution determined by this method has high accuracy.

Step 2: Extract constant value parameters as the state quantities of the system. Considering that the actual observed value η is the quaternion of the observed coordinate system relative to the inertial system, the relationship between it and λ

η=λομ

Among them, "ο" represents the quaternion multiplication operation, and μ is the attitude quaternion between the observation coordinate system and the main axis coordinate system of the body. Since the quaternion multiplication does not involve high-order term operations, and μ is a constant, each item of η can be a linear combination of each item of λ. From the form of formula (3), we can see that the observed quantity η also has periodic properties, that is,

η=B _4×6 ξ (5)

Where B is a constant matrix satisfying Aξ=Bξομ. Expand the matrix B by row to get the row vector b _24×1 , and consider the constant parameter to be identified and θ, together constitute the constant state quantity x of the system

b, θ, are constant parameters to be identified;

Step 3: Construct an extended Kalman filter.

The observation equation of the system is

h(x)=η

Calculate the partial derivative of h relative to the state quantity x to obtain the observation sensitivity matrix

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

Φ＝Ι _14×14 (7)

Wherein, I _14×14 is an identity matrix of 14×14;

The error of the state quantity caused by the disturbance torque is considered to be white noise of Gaussian distribution, then the variance matrix Q _k of the system process noise is defined as

Q _k = σ ² Ι _14×14 (8)

Among them, ^σ2 is the distribution variance of the process noise;

Substituting equations (6)-(8) into the general equation of the Kalman filter is the step of identifying the parameters of the asymmetric tumbling target. Similar to the case of the symmetric target, the items of the b vector in the state quantity satisfy the linear form, and there is no need to give an accurate initial value. while θ and Involves trigonometric function operations, so it is necessary to use discrete Fourier transform to determine the characteristic angular frequency θ and initial value.

Step 4: After the estimated value of the state quantity x is obtained through the Kalman filter, the estimated value is substituted into formula (4) and the least square solution is obtained, and the matrix D and the estimation of the angular velocity parameters ω _xm , ω _ym and ω _zm can be obtained value. Further, the estimated value of the moment of inertia parameter is calculated by the following formula

p _x , p _y , p _z are the moment of inertia ratio parameters to be identified, and K is the first kind of complete elliptic integral coefficient corresponding to the elliptic function of angular velocity.

In order to better illustrate the purpose and advantages of the present invention, the contents of the present invention will be further described below in conjunction with the accompanying drawings and examples:

Assuming that the tumbling target is floating freely in space, the function of the four variables of the attitude quaternion changing with time can be measured by using a stereo vision device or a laser range finder, as shown in Figure 2. Due to the influence of disturbance torque and observation error, the measurement result is polluted by noise. By applying this method, the inertial parameters p _x , p _y and p _z of the tumbling satellite can be estimated in real time by using these observation noises, which specifically includes the following steps:

Step 1: first use the fast Fourier transform algorithm to transform part of the data of λ ₀ (t) into the frequency domain, as shown in FIG. 3 . There are three peaks in the frequency domain, assign the value of the abscissa of the middle peak to θ, and assign the interval between the abscissas of the peaks to as the initial value. The initial values of other parameters in x are assigned 0.

Step 2: Taking x as the state parameter and taking the real-time observation of the attitude quaternion λ as input, construct a Kalman filter to gradually estimate more accurate values of the state parameter. As shown in Figure 4, the errors of the attitude quaternion and angular velocity calculated and predicted by the estimated value of the constant parameter are smaller than the error of the direct observation value.

Step 3: Calculate the values of the rolling satellite inertia parameters p _x , p _y and p _z by using the values of each parameter x. The curves of the estimated values changing with time are shown in Figure 5. It can be seen that the estimated values of the inertial parameters increase with the increase of observations approaching the real value.

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

Table 1. System parameter values used in the example

Claims (5)

1. A non-contact inertial coefficient identification method for tumbling satellites with arbitrary shapes, characterized in that it comprises the steps: 1) establish the approximate solution of the motion equation of the rolling target quaternion λ by an approximate method; 2) Extract the constant parameter as the state quantity x of the system through the approximate solution; 3) Construct an extended Kalman filter; 4) Obtain the estimated value of the state quantity of the system through the Kalman filter, and use the estimated value to estimate the inertial parameters of the rolling satellite. 2. The non-contact inertia coefficient identification method of a tumbling satellite of any shape according to claim 1, wherein step 1) is specifically: 1-1) Establish the kinematics model of the attitude quaternion of the rolling target: $\begin{array}{c}{\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; y</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; z</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; )</span>\\ {\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; z</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; y</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; )</span>\\ {\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; y</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; z</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; )</span>\\ {\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; z</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; the\; y</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>\end{array}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$ Among them, λ ₀ , λ ₁ , λ ₂ , λ ₃ are the four components of the attitude quaternion of the main axis coordinate system of the rolling target relative to the inertial coordinate system, and ω _x , ω _y , ω _z are the angular velocity of the target in the body coordinates projections in the system; 1-2) Neglecting the small amount of high frequency, the approximate solution of the angular velocity of the rolling target is obtained: Among them, ω _xm , ω _ym , ω _zm are the maximum values of the angular velocity components in the directions of the three coordinate axes of the principal axis coordinate system of the body respectively, t is the time, is one of the constant parameters to be identified, and ε _x , ε _y , ε _z are all high-order infinitesimal quantities in the direction of the coordinate axis; 1-3) Use the sine function containing three characteristic frequencies to construct the approximate solution of the quaternion differential equation describing the attitude motion, the specific form of the solution of the quaternion λ=[λ ₀ λ ₁ λ ₂ λ ₃ ] ^T is: λ＝A _4×6 ξ (3) where A _4×6 is a constant matrix with 4 rows and 6 columns, and Among them, θ is also a constant value parameter to be identified; 1-4) Transform the differential equation of the system into an algebraic equation, substitute formula (2) and formula (3) into formula (1) at the same time, and use the trigonometric function operation rules to combine the same frequency items. From the regularity of the trigonometric function, if Equation (1) is established, and the coefficients before each frequency must be zero, then an equation system consisting of 56 arithmetic equations is obtained: Where a is the row vector generated after the matrix A _4×6 is expanded into a vector by row, for the reason The constant coefficient matrix determined by θ and ω _xm , ω _ym , ω _zm , the specific form of each item is derived from formula (3), and the subscript indicates the dimension of the matrix. Considering the initial value λ(0) of the system, the upper expands to: Among them, the value of the constant matrix D' is deduced by D'a=λ(0): ${\mathrm{<span\; class="notranslate">\; D.</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; =</span>\left[\begin{array}{cccccccccccccccccccccccc}\mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& & & & & & & & & & & & & & & & & & \\ & & & & & & \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& & & & & & & & & & & & \\ & & & & & & & & & & & & \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& & & & & & \\ & & & & & & & & & & & & & & & & & & \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right]$ Equation (4) is an overdetermined equation system with the number of constraints greater than the number of variables, using the least squares criterion to obtain a, And the least squares solution of θ, the approximate solution of the quaternion λ about time t is obtained, that is, the approximate equation of the attitude quaternion changing with time, that is, the approximate equation of the attitude motion of the rolling target. 3. The non-contact method for identifying the coefficient of inertia of a tumbling satellite of any shape according to claim 2, wherein step 2) is specifically: The actual observation value η is the quaternion of the observation coordinate system relative to the inertial system, and the relationship between it and λ is: in Indicates the quaternion multiplication operation; μ is the attitude quaternion between the observation coordinate system and the main axis coordinate system of the body; Since the quaternion multiplication does not involve high-order term operations, and μ is a constant, the items of η are linear combinations of the items of λ. It can be seen from formula (3) that the observed quantity η also has periodic properties, that is: η=B _4×6 ξ (5) where B is satisfied The constant matrix of , expand the matrix B by row to get the row vector b _24×1 , and consider the constant parameter to be identified and θ, together constitute the constant state quantity of the system: in, are constant parameters to be identified. 4. The non-contact method for identifying the coefficient of inertia of a tumbling satellite of any shape according to claim 2, wherein step 3) is specifically: The observation equation of the system is h(x)=η Calculate the partial derivative of h relative to the state quantity x to obtain the observation sensitivity matrix: Since the state quantities of the system are constant parameters, the state transition matrix Φ is: Φ＝Ι _14×14 (7) Wherein, I _14×14 is the identity matrix of 14×14; The error of the state quantity caused by the disturbance torque is considered to be white noise of Gaussian distribution, then the variance matrix Q _k of the system process noise is defined as: Q _k = σ ² Ι _14×14 (8) Among them, ^σ2 is the distribution variance of the process noise; Substituting equations (6)-(8) into the general equation of the Kalman filter is the step of identifying the parameters of the asymmetric tumbling target, and constructing a generalized Kalman filter composed of the parameters to be identified as state variables. 5. The non-contact method for identifying the coefficient of inertia of a tumbling satellite of any shape according to claim 4, wherein step 4) is specifically: After the estimated value of the state quantity x is obtained through the Kalman filter, the estimated value is substituted into formula (4) and the least square solution is obtained, that is, the estimated value of the matrix D and the angular velocity parameters ω _xm , ω _ym and ω _zm is obtained, and further Solving for an estimate of the inertia parameter: Among them, p _x , p _y , p _z are the moment of inertia ratio parameters to be identified, and K is the first kind of complete elliptic integral coefficient corresponding to the elliptic function of angular velocity.