CN103940433A - Satellite attitude determining method based on improved self-adaptive square root UKF (Unscented Kalman Filter) algorithm - Google Patents

Abstract

The invention relates to a satellite attitude determining method based on an improved self-adaptive square root UKF (Unscented Kalman Filter) algorithm, which belongs to the technical field of satellite attitude determination and solves the problems that a satellite attitude determining system is unstable, the satellite attitude precision is low, and the traceability to an virtual condition of a satellite is weak due to overlarge round-off errors numerically calculated by virtue of existing EKF (Extended Kalman Filter), UKF (Unscented Kalman Filter) and SRUKF (Square Root Unscented Kalman Filter) algorithms when the satellite attitude determining system suffers from uncertain interferences and is influenced by noise. The satellite attitude determining method comprises the main realization processes: estimating an error quaternion and a gyroscopic drift error by virtue of the improved self-adaptive square root UKF; substituting gyroscopic measurement value and the estimated gyroscopic drift error into an attitude kinematical equation to calculate an attitude quaternion; correcting the calculated attitude quaternion by virtue of the estimated error quaternion; carrying out attitude resolving by virtue of the corrected attitude quaternion so as to determine the attitude of the satellite. The satellite attitude determining method is suitable for the technical field of the satellite attitude determination.

Description

A kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm

Technical field

The present invention relates to the high-precision satellite attitude determination method of star sensor and gyro, belong to Satellite attitude determination technical field.

Background technology

The attitude measurement system of star sensor and gyro composition, because its precision is higher, thereby is widely used in Satellite Attitude Determination System.For the nonlinear system of its composition, adopt nonlinear filtering technique to determine that the method for attitude is widely used.Extended Kalman filter (EKF) is because its method is simple, and the advantages such as easy realization are widely used in the engineering problem of processing nonlinear estimation.But EKF does first approximation processing to nonlinear equation, ignores all the other higher order terms, thereby nonlinear problem is converted into linear problem.Non-linear when stronger when system, EKF runs counter to local linear hypothesis, and uncared-for higher order term can bring large error, causes EKF arithmetic accuracy to decline, and even causes the filtering divergence of EKF algorithm.EKF need to calculate Jacobi (Jacobian) matrix in the time of linearization process in addition, its computation process very complicated and easily makeing mistakes.

For the problems referred to above, the people such as Julier have proposed Wuji Kalman filter (UKF) algorithm.With respect to EKF, UKF adopts UT transfer pair nonlinear probability Density Distribution to be similar to, and has and does not need to calculate Jacobian matrix, and estimated accuracy is higher, obtains widespread use in recent years in the middle of attitude of flight vehicle estimation problem.But UKF tends to exist round-off error in numerical evaluation, may destroy nonnegative definiteness and the symmetry of system estimation error covariance matrix, cause convergence of algorithm speed slow, even cause the unstable of algorithm.Can solve preferably for this problem square root UKF (SRUKF) filtering algorithm, the method is used for reference square root in kalman filtering and is decomposed filter thought, in the process of filtering, adopt the QR of matrix to decompose and the square root of covariance matrix is directly propagated and upgraded to Cholesky decomposition result, solve the negative definiteness problem of the error covariance that in UKF algorithm, the error of calculation may cause, improved counting yield and the numerical stability of filtering algorithm.But, when numerical value cumulative errors is too large or weights w ₀ ^c selecting inappropriate time, still there is filtering divergence problem in the square root UKF of standard.In addition, EKF, UKF and SRUKF all require system model accurately and noise statistics known.In the time that system exists uncertain interference and noise contributions, said method does not possess good estimated accuracy, tracking power and robustness.

Summary of the invention

The object of the invention is to propose a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm, to solve in the time that Satellite Attitude Determination System is subject to affecting of uncertain interference and noise, due to existing EKF, the too large caused Satellite Attitude Determination System of round-off error of UKF and SRUKF algorithm numerical evaluation unstable and the problem such as a little less than low and satellite virtual condition tracking power to the precision of Satellite Attitude Estimation.

The present invention for solving the problems of the technologies described above adopted technical scheme is:

A kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm of the present invention, realizes according to following steps: step 1, set up gyro to measure model; Step 2, set up satellite attitude kinematics equation; The system state equation of step 3, the state variable of foundation based on error quaternion and gyroscopic drift error composition; Step 4, set up error system observation equation; Step 5, utilize improved self-adaptation square root UKF evaluated error hypercomplex number and gyroscopic drift error; Step 6, the gyroscopic drift error substitution attitude motion equation calculating attitude quaternion that utilizes gyro to measure value and estimate; The error quaternion that step 7, utilization estimate is revised the attitude quaternion calculating; The attitude quaternion that step 8, utilization are revised carries out attitude algorithm, determines the attitude of satellite.

The invention has the beneficial effects as follows:

One, the present invention has improved the stability due to unstable the caused Satellite Attitude Determination System of existing algorithm, makes the stability of Satellite Attitude Determination System better.

Two, on improved square root UKF basis, introduce adaptive factor μ _k+1, make square root UKF there is adaptivity, thereby in the situation that system has unknown disturbances, improved Satellite Attitude Estimation precision.

Three, the present invention has improved UKF and the tracking power of square root UKF to unknown disturbances and mutation status, satellite virtual condition tracking power is strengthened, thereby also can estimate preferably the attitude of satellite when uncertain at model uncertainty and noise statistics, therefore more applicable in the middle of environment complicated and changeable.

Four, with respect to STF method, the present invention (IASRUKF) is improving 66%～67% aspect the estimated accuracy of roll angle, aspect the estimated accuracy of the angle of pitch, improving 68%～69%, aspect the estimated accuracy of crab angle, improving 58%～%59; With respect to SRUKF method, the present invention is improving 97%～98% aspect the estimated accuracy of roll angle, aspect the estimated accuracy of the angle of pitch, is improving 97%～98%, aspect the estimated accuracy of crab angle, is improving 97%～98%.

Brief description of the drawings

Fig. 1 is process flow diagram of the present invention; Fig. 2 is the improved self-adaptation square root UKF algorithm flow chart of attitude of satellite error system; Fig. 3 is the roll angle graph of errors comparison diagram of STF (strong tracking filter), SRUKF (square root UKF) and IASRUKF of the present invention; Fig. 4 is the angle of pitch graph of errors comparison diagram of STF (strong tracking filter), SRUKF (square root UKF) and IASRUKF of the present invention; Fig. 5 is the crab angle graph of errors comparison diagram of STF (strong tracking filter), SRUKF (square root UKF) and IASRUKF of the present invention.

Embodiment

Embodiment one: a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm described in present embodiment, is characterized in that said method comprising the steps of:

Step 1, set up gyro to measure model;

Step 2, set up satellite attitude kinematics equation;

The system state equation of step 3, the state variable of foundation based on error quaternion and gyroscopic drift error composition;

Step 4, set up error system observation equation;

Step 5, utilize improved self-adaptation square root UKF evaluated error hypercomplex number and gyroscopic drift error;

Step 6, the gyroscopic drift error substitution attitude motion equation calculating attitude quaternion that utilizes gyro to measure value and estimate;

The error quaternion that step 7, utilization estimate is revised the attitude quaternion calculating;

The attitude quaternion that step 8, utilization are revised carries out attitude algorithm, determines the attitude of satellite.

Embodiment two: present embodiment is different from embodiment one: the detailed process of setting up gyro to measure model described in step 1 is: be the same coordinate system at the coordinate of gyro to measure coordinate system and celestial body, gyro to measure model is

$\left\{\begin{array}{c}g\left(t\right)=\mathrm{\ω}\left(t\right)+\mathrm{\β}\left(t\right)+{\mathrm{\η}}_u\left(t\right)\\ \stackrel{\·}{\mathrm{\β}\left(t\right)}={\mathrm{\η}}_v\left(t\right)\end{array}\right.---\left(1\right)$

In formula, the measurement output valve that g (t) is gyro, ω (t) is gyro true angular velocity, β (t) is gyroscopic drift, η _uand η (t) _v(t) be mutual incoherent Gaussian noise, meet:

$\left\{\begin{array}{c}E\left\{{\mathrm{\η}}_u\left(t\right){\mathrm{\η}}_u^T\left(\mathrm{\τ}\right)\right\}={\mathrm{\σ}}_u^2\mathrm{\δ}(t-\mathrm{\τ})I_{3\×3}\\ E\left\{{\mathrm{\η}}_v\left(t\right){\mathrm{\η}}_v^T\left(\mathrm{\τ}\right)\right\}={\mathrm{\σ}}_v^2\mathrm{\δ}(t-\mathrm{\τ})I_{3\×3}\end{array}\right.---\left(2\right)$

In formula, with be respectively white noise, η _uand η (t) _v(t) be mean square deviation, δ (t-τ) is Dirac function.Other step and parameter are identical with embodiment one.

Embodiment three: present embodiment is different from embodiment one or two: the detailed process of setting up satellite attitude kinematics equation described in step 2 is: attitude of satellite hypercomplex number is defined as

q＝[q ₁ q ₂ q ₃ q ₄] ^T (3)

In formula, q ₄=cos (θ/2), with the θ unit of being respectively rotating vector and rotation angle;

Hypercomplex number meets following constraint:

$q_1^2+q_2^2+q_3^2+q_4^2=1---\left(4\right)$

Represent that by hypercomplex number satellite attitude kinematics equation is:

$\frac{\mathrm{dq}\left(t\right)}{\mathrm{dt}}=\frac{1}{2}\mathrm{\Ω}\left(\mathrm{\ω}\left(t\right)\right)q\left(t\right)---\left(5\right)$

In formula, the true angular velocity that ω (t) is satellite, and

$\mathrm{\Ω}\left(\mathrm{\ω}\left(t\right)\right)={\left[\begin{array}{cc}-{[\mathrm{\ω}\left(t\right)\×]}_{3\×3}& \mathrm{\ω}\left(t\right)\\ -{\mathrm{\ω}}^T\left(t\right)& 0\end{array}\right]}_{4\×4},[\mathrm{\ω}\left(t\right)\×]=\left[\begin{array}{ccc}0& -{\mathrm{\ω}}_z\left(t\right)& {\mathrm{\ω}}_y\left(t\right)\\ {\mathrm{\ω}}_z\left(t\right)& 0& -{\mathrm{\ω}}_x\left(t\right)\\ -{\mathrm{\ω}}_y\left(t\right)& {\mathrm{\ω}}_x\left(t\right)& 0\end{array}\right].$ Other step and parameter are identical with embodiment one or two.

Embodiment four: present embodiment is different from one of embodiment one to three: the detailed process of the system state equation of the state variable of the foundation described in step 3 based on error quaternion and gyroscopic drift error composition is: hypercomplex number exists norm constraint, if selected using four quaternion components as state variable, variance matrix is unusual.Therefore for fear of the unusual problem of variance matrix, adopt multiplicative quaternion to define true hypercomplex number and hypercomplex number calculated value between error quaternion be state variable:

$\mathrm{\δq}={\hat{q}}^{-1}\⊗q=\left[\begin{array}{cc}\mathrm{\δe}& \mathrm{\δq}\end{array}\right]---\left(6\right)$

Get the vector part δ e of error quaternion and gyroscopic drift error △ b as state variable, the state variable of system is △ x=[δ e △ b], had by satellite attitude kinematics equation:

$\stackrel{\·}{q}=\frac{1}{2}q\⊗\mathrm{\ω}---\left(7\right)$

$\stackrel{\stackrel{\·}{^}}{q}=\frac{1}{2}\hat{q}\⊗\widehat{\mathrm{\ω}}---\left(8\right)$

In formula, it is the calculated value of Satellite Angle speed;

To formula (6) differentiate, and by formula (7) and formula (8) substitution formula (6):

$\mathrm{\δ}\stackrel{\·}{q}=\frac{1}{2}[\mathrm{\δq}\⊗\mathrm{\ω}-\widehat{\mathrm{\ω}}\⊗\mathrm{\δq}]---\left(9\right)$

Definition:

$\mathrm{\δ\ω}=[\mathrm{\ω}-\widehat{\mathrm{\ω}}]=[(b-\hat{b})+{\mathrm{\η}}_u]---\left(10\right)$

Formula (10) substitution formula (9) is obtained:

$\mathrm{\δ}\stackrel{\·}{q}=\frac{1}{2}[\mathrm{\δq}\⊗\widehat{\mathrm{\ω}}-\widehat{\mathrm{\ω}}\⊗\mathrm{\δq}]+\frac{1}{2}\mathrm{\δq}\⊗\mathrm{\δ\ω}---\left(11\right)$

By formula (11) linearization, have:

$\mathrm{\δ}\stackrel{\·}{q}=-\widehat{\mathrm{\ω}}\×\mathrm{\δe}-\frac{1}{2}(\mathrm{\Δb}+{\mathrm{\η}}_u)---\left(12\right)$

The system state equation that can obtain error state variable is:

$\mathrm{\Δ}\stackrel{\·}{X}\left(t\right)=F\left(t\right)\mathrm{\ΔX}\left(t\right)+G\left(t\right)\mathrm{\ϵ}\left(t\right)---\left(13\right)$

In formula,

$F\left(t\right)=\left[\begin{array}{cc}-[\mathrm{\ω}\left(t\right)\×]& \frac{1}{2}{I}_{3\×3}\\ 0_{3\×3}& 0_{3\×3}\end{array}\right],G\left(t\right)=\left[\begin{array}{cc}-\frac{1}{2}{I}_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& I_{3\×3}\end{array}\right],\mathrm{\ϵ}\left(t\right)=\begin{array}{c}\left[\begin{array}{c}{\mathrm{\η}}_u\\ 0\end{array}\right]\end{array}.$ Other step and parameter are identical with one of embodiment one to three.

Embodiment five: present embodiment is different from one of embodiment one to four: the detailed process of setting up error system observation equation described in step 4 is: be using hypercomplex number q as output quantity due to what adopt, the hypercomplex number error that star sensor is obtained is as measuring output, get its vector part as observed quantity, have observation equation:

z＝Hx+v (14)

In formula, H is observing matrix, and v is observation noise;

$H=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right]$

Adopt fourth order Runge-Kutta method to carry out discretize to the attitude error system equation of formula (13) and formula (14) composition.Other step and parameter are identical with one of embodiment one to four.

Embodiment six: present embodiment is different from one of embodiment one to five: the detailed process of utilizing improved self-adaptation square root UKF evaluated error hypercomplex number and gyroscopic drift error described in step 5 is:

(1) init state state error covariance matrix square root S ₀for

${\hat{X}}_0=E\left[X_0\right]---\left(15\right)$

$S_0=\mathrm{chol}\left(E\right[(X_0-{\hat{X}}_0{(X_0-{\hat{X}}_0)}^T)\left]\right)---\left(16\right)$

(2) for k=1,2,3 ..., n performing step is as follows:

1. calculate sigma point

${\mathrm{\ξ}}_k=\left[\begin{array}{ccc}{\hat{X}}_k& {\hat{X}}_k+\mathrm{\γ}{S}_k& {\hat{X}}_k-\mathrm{\γ}{S}_k\end{array}\right]---\left(17\right)$

In formula, λ=α ²(n+ κ)-n, the value of α is between 0.001～1, and n is the dimension of system state vector, and κ is the 3rd scale factor, and value is 0;

2. the time upgrades

The nonlinear propagation of Sigma sampled point:

ξ _i,k+1|k＝F _k+1|kξ _i,k (18)

The one-step prediction value of state value is calculated:

${\hat{X}}_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^m{\mathrm{\ξ}}_{k+1|k}---\left(19\right)$

The one-step prediction of observed reading calculates:

${\hat{Z}}_{k+1|k}=H{\hat{X}}_{k+1|k}---\left(20\right)$

The residual computations of predicted value and actual value:

$e_{k+1}=Z_{k+1}-{\hat{Z}}_{k+1|k}---\left(21\right)$

3. calculate self-adaptation fading factor μ _k+1

${\mathrm{\&mu;}}_{k+1,i}=\left\{\begin{array}{c}{\mathrm{\&mu;}}_{k+1,i},{\mathrm{\&mu;}}_{k+1,i}\&GreaterEqual;1\\ 1,{\mathrm{\&mu;}}_{k+1,i}<1\end{array}\right.,i=\mathrm{1,2},...,n---\left(22\right)$

μ _k+1,i＝trace(N _k+1)/trace(M _ii,k+！) (23)

N _k+1＝V _k+1-HQ _k+1H ^T-lR _k+1 (24)

M _k+1＝J _k+1HH ^T (25)

$J_{k+1}=\underset{i=1}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^c({\mathrm{\&xi;}}_{k+1|k,i}-{\mathrm{\&xi;}}_{0,k+1|k}){({\mathrm{\&xi;}}_{k+1|k,i}-{\mathrm{\&xi;}}_{0,k+1|k})}^T---\left(26\right)$

$V_{k+1}=\left\{\begin{array}{c}{e}_{k+1}{e}_{k+1}^T,k=0\\ \frac{\mathrm{\&rho;}{V}_k+e_{k+1}{e}_{k+1}^T}{1+\mathrm{\&rho;}},k\&GreaterEqual;1\end{array}\right.---\left(27\right)$

4. measure and upgrade

The square root of state one-step prediction covariance matrix calculates:

$S_{k+1|k}=\mathrm{qr}\left\{\right[\mathrm{diag}\left({\mathrm{\&mu;}}_{k+1}\right)\sqrt{W_1^c}({\mathrm{\&xi;}}_{1:2n,k+1|k}-{\mathrm{\&xi;}}_{0,k+1|k}),G_k\sqrt{Q_k}\left]\right\}---\left(28\right)$

Cross-covariance calculates:

$P_{\mathrm{XZ},k+1}=\left(S_{k+1|k}{S}_{k+1|k}^T\right)H^T---\left(29\right)$

The square root of output covariance matrix calculates:

U _k+1＝HS _k+1|k (30)

$S_{Z,k+1}=\mathrm{qr}\left\{\right[U_{k+1},\sqrt{R_{k+1}}\left]\right\}---\left(31\right)$

State gain matrix calculates:

$K_{k+1}=(P_{\mathrm{XZ},k+1}/S_{Z,k+1}^T)/S_{Z,k+1}---\left(32\right)$

State estimation value is calculated:

${\hat{X}}_{k+1|k+1}={\hat{X}}_{k+1|k}+K_{k+1}{e}_{k+1}---\left(33\right)$

The square root of state error covariance matrix calculates:

$S_{k+1}=S_{k+1|k}-S_{k+1|k}{U}_{k+1}^T{S}_{z,k+1}^{-T}{(S_{z,k+1}+\sqrt{R_{k+1}})}^{-1}{U}_{k+1}---\left(34\right)$

The calculation of parameter of using in said process is as follows:

$W_0^m=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}---\left(35\right)$

$W_0^c=W_0^m+(1-{\mathrm{\&alpha;}}^2+\mathrm{\&beta;})---\left(36\right)$

other step and parameter are identical with one of embodiment one to five.

Embodiment seven: present embodiment is different from one of embodiment one to six: the attitude algorithm method described in step 8 is:

θ＝arcsin(-2(q ₁q ₃-q ₂q ₄)) (39)

$\mathrm{\&psi;}=\arctan (2(q_2{q}_1+q_3{q}_4)/(1-2(q_3^2+q_2^2))---\left(40\right)$

In formula, θ, ψ represents respectively roll angle, the angle of pitch and crab angle.Other step and parameter are identical with one of embodiment one to six.

For advantage of the present invention is described, the present invention to be compared with strong tracking filter (STF) method and square root UKF (SRUKF) method, simulation parameter is set to: the initial attitude angular velocity of satellite is: ω ₀=10 ^-2× [2 2 2] ^tdeg/s, gyroscopic drift is: β ₀=[0.05 0.05 0.05] ^tdeg/s, the measurement noise of gyro and the mean square deviation of drift noise are respectively: σ _u=0.5 (°)/h and σ _v=0.04 (°)/h, the mean square deviation of measuring noise is: σ _q=10 ", the filtering sampling time is △ t=1s, when simulation time is 200s, adds mutation disturbance: △ x=[0,0,0,0.02,0.02,0.02] ^t, when simulation time is 500s, the noise statistics of system is expanded to 20 times.Obtain the simulation experiment result as shown in Fig. 3～Fig. 5.By the analysis to the simulation experiment result, estimated accuracy of the present invention is best, with respect to STF method, the present invention (IASRUKF) is improving 66%～67% aspect the estimated accuracy of roll angle, aspect the estimated accuracy of the angle of pitch, improving 68%～69%, aspect the estimated accuracy of crab angle, improving 58%～%59; With respect to SRUKF method, the present invention is improving 97%～98% aspect the estimated accuracy of roll angle, aspect the estimated accuracy of the angle of pitch, is improving 97%～98%, aspect the estimated accuracy of crab angle, is improving 97%～98%.

Claims (7)

1. the Satellite attitude determination method based on improved self-adaptation square root UKF algorithm, is characterized in that said method comprising the steps of:

Step 1, set up gyro to measure model;

Step 2, set up satellite attitude kinematics equation;

The system state equation of step 3, the state variable of foundation based on error quaternion and gyroscopic drift error composition;

Step 4, set up error system observation equation;

Step 5, utilize improved self-adaptation square root UKF evaluated error hypercomplex number and gyroscopic drift error;

Step 6, the gyroscopic drift error substitution attitude motion equation calculating attitude quaternion that utilizes gyro to measure value and estimate;

The error quaternion that step 7, utilization estimate is revised the attitude quaternion calculating;

The attitude quaternion that step 8, utilization are revised carries out attitude algorithm, determines the attitude of satellite.

2. a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm according to claim 1, it is characterized in that the detailed process of setting up gyro to measure model described in step 1 is: be the same coordinate system at the coordinate of gyro to measure coordinate system and celestial body, gyro to measure model is

$\left\{\begin{array}{c}g\left(t\right)=\mathrm{\&omega;}\left(t\right)+\mathrm{\&beta;}\left(t\right)+{\mathrm{\&eta;}}_u\left(t\right)\\ \stackrel{\&CenterDot;}{\mathrm{\&beta;}\left(t\right)}={\mathrm{\&eta;}}_v\left(t\right)\end{array}\right.---\left(1\right)$

In formula, the measurement output valve that g (t) is gyro, ω (t) is gyro true angular velocity, β (t) is gyroscopic drift, η _uand η (t) _v(t) be mutual incoherent Gaussian noise, meet:

$\left\{\begin{array}{c}E\left\{{\mathrm{\&eta;}}_u\left(t\right){\mathrm{\&eta;}}_u^T\left(\mathrm{\&tau;}\right)\right\}={\mathrm{\&sigma;}}_u^2\mathrm{\&delta;}(t-\mathrm{\&tau;})I_{3\&times;3}\\ E\left\{{\mathrm{\&eta;}}_v\left(t\right){\mathrm{\&eta;}}_v^T\left(\mathrm{\&tau;}\right)\right\}={\mathrm{\&sigma;}}_v^2\mathrm{\&delta;}(t-\mathrm{\&tau;})I_{3\&times;3}\end{array}\right.---\left(2\right)$

In formula, with be respectively white noise, η _uand η (t) _v(t) be mean square deviation, δ (t-τ) is Dirac function.

3. a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm according to claim 2, is characterized in that the detailed process of setting up satellite attitude kinematics equation described in step 2 is: attitude of satellite hypercomplex number is defined as

q＝[q ₁ q ₂ q ₃ q ₄] ^T (3)

In formula, q ₄=cos (θ/2), with the θ unit of being respectively rotating vector and rotation angle;

Hypercomplex number meets following constraint:

$q_1^2+q_2^2+q_3^2+q_4^2=1---\left(4\right)$

Represent that by hypercomplex number satellite attitude kinematics equation is:

$\frac{\mathrm{dq}\left(t\right)}{\mathrm{dt}}=\frac{1}{2}\mathrm{\&Omega;}\left(\mathrm{\&omega;}\left(t\right)\right)q\left(t\right)---\left(5\right)$

In formula, the true angular velocity that ω (t) is satellite, and

$\mathrm{\&Omega;}\left(\mathrm{\&omega;}\left(t\right)\right)={\left[\begin{array}{cc}-{[\mathrm{\&omega;}\left(t\right)\&times;]}_{3\&times;3}& \mathrm{\&omega;}\left(t\right)\\ -{\mathrm{\&omega;}}^T\left(t\right)& 0\end{array}\right]}_{4\&times;4},[\mathrm{\&omega;}\left(t\right)\&times;]=\left[\begin{array}{ccc}0& -{\mathrm{\&omega;}}_z\left(t\right)& {\mathrm{\&omega;}}_y\left(t\right)\\ {\mathrm{\&omega;}}_z\left(t\right)& 0& -{\mathrm{\&omega;}}_x\left(t\right)\\ -{\mathrm{\&omega;}}_y\left(t\right)& {\mathrm{\&omega;}}_x\left(t\right)& 0\end{array}\right].$

4. a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm according to claim 3, is characterized in that the detailed process of the system state equation of foundation described in the step 3 state variable based on error quaternion and gyroscopic drift error composition is: adopt multiplicative quaternion to define true hypercomplex number and hypercomplex number calculated value between error quaternion be state variable:

$\mathrm{\&delta;q}={\hat{q}}^{-1}\&CircleTimes;q=\left[\begin{array}{cc}\mathrm{\&delta;e}& \mathrm{\&delta;q}\end{array}\right]---\left(6\right)$

Get the vector part δ e of error quaternion and gyroscopic drift error △ b as state variable, the state variable of system is △ x=[δ e △ b], had by satellite attitude kinematics equation:

$\stackrel{\&CenterDot;}{q}=\frac{1}{2}q\&CircleTimes;\mathrm{\&omega;}---\left(7\right)$

$\stackrel{\stackrel{\&CenterDot;}{^}}{q}=\frac{1}{2}\hat{q}\&CircleTimes;\widehat{\mathrm{\&omega;}}---\left(8\right)$

In formula, it is the calculated value of Satellite Angle speed;

To formula (6) differentiate, and by formula (7) and formula (8) substitution formula (6):

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{q}=\frac{1}{2}[\mathrm{\&delta;q}\&CircleTimes;\mathrm{\&omega;}-\widehat{\mathrm{\&omega;}}\&CircleTimes;\mathrm{\&delta;q}]---\left(9\right)$

Definition:

$\mathrm{\&delta;\&omega;}=[\mathrm{\&omega;}-\widehat{\mathrm{\&omega;}}]=[(b-\hat{b})+{\mathrm{\&eta;}}_u]---\left(10\right)$

Formula (10) substitution formula (9) is obtained:

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{q}=\frac{1}{2}[\mathrm{\&delta;q}\&CircleTimes;\widehat{\mathrm{\&omega;}}-\widehat{\mathrm{\&omega;}}\&CircleTimes;\mathrm{\&delta;q}]+\frac{1}{2}\mathrm{\&delta;q}\&CircleTimes;\mathrm{\&delta;\&omega;}---\left(11\right)$

By formula (11) linearization, have:

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{q}=-\widehat{\mathrm{\&omega;}}\&times;\mathrm{\&delta;e}-\frac{1}{2}(\mathrm{\&Delta;b}+{\mathrm{\&eta;}}_u)---\left(12\right)$

The system state equation that can obtain error state variable is:

$\mathrm{\&Delta;}\stackrel{\&CenterDot;}{X}\left(t\right)=F\left(t\right)\mathrm{\&Delta;X}\left(t\right)+G\left(t\right)\mathrm{\&epsiv;}\left(t\right)---\left(13\right)$

In formula,

$F\left(t\right)=\left[\begin{array}{cc}-[\mathrm{\&omega;}\left(t\right)\&times;]& \frac{1}{2}{I}_{3\&times;3}\\ 0_{3\&times;3}& 0_{3\&times;3}\end{array}\right],G\left(t\right)=\left[\begin{array}{cc}-\frac{1}{2}{I}_{3\&times;3}& 0_{3\&times;3}\\ 0_{3\&times;3}& I_{3\&times;3}\end{array}\right],\mathrm{\&epsiv;}\left(t\right)=\begin{array}{c}\left[\begin{array}{c}{\mathrm{\&eta;}}_u\\ 0\end{array}\right]\end{array}.$

5. a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm according to claim 4, it is characterized in that the detailed process of setting up error system observation equation described in step 4 is: is using hypercomplex number q as output quantity due to what adopt, the hypercomplex number error that star sensor is obtained is as measuring output, get its vector part as observed quantity, have observation equation:

z＝Hx+v (14)

In formula, H is observing matrix, and v is observation noise;

$H=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right]$

Adopt fourth order Runge-Kutta method to carry out discretize to the attitude error system equation of formula (13) and formula (14) composition.

6. a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm according to claim 5, is characterized in that the detailed process of utilizing improved self-adaptation square root UKF evaluated error hypercomplex number and gyroscopic drift error described in step 5 is:

(1) init state state error covariance matrix square root S ₀for

${\hat{X}}_0=E\left[X_0\right]---\left(15\right)$

$S_0=\mathrm{chol}\left(E\right[(X_0-{\hat{X}}_0{(X_0-{\hat{X}}_0)}^T)\left]\right)---\left(16\right)$

(2) for k=1,2,3 ..., n performing step is as follows:

1. calculate sigma point

${\mathrm{\&xi;}}_k=\left[\begin{array}{ccc}{\hat{X}}_k& {\hat{X}}_k+\mathrm{\&gamma;}{S}_k& {\hat{X}}_k-\mathrm{\&gamma;}{S}_k\end{array}\right]---\left(17\right)$

In formula, λ=α ²(n+ κ)-n, the value of α is between 0.001～1, and n is the dimension of system state vector, and κ is the 3rd scale factor, and value is 0;

2. the time upgrades

The nonlinear propagation of Sigma sampled point:

ξ _i,k+1|k＝F _k+1|kξ _i,k (18)

The one-step prediction value of state value is calculated:

${\hat{X}}_{k+1|k}=\underset{i=0}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^m{\mathrm{\&xi;}}_{k+1|k}---\left(19\right)$

The one-step prediction of observed reading calculates:

${\hat{Z}}_{k+1|k}=H{\hat{X}}_{k+1|k}---\left(20\right)$

The residual computations of predicted value and actual value:

$e_{k+1}=Z_{k+1}-{\hat{Z}}_{k+1|k}---\left(21\right)$

3. calculate self-adaptation fading factor μ _k+1

${\mathrm{\&mu;}}_{k+1,i}=\left\{\begin{array}{c}{\mathrm{\&mu;}}_{k+1,i},{\mathrm{\&mu;}}_{k+1,i}\&GreaterEqual;1\\ 1,{\mathrm{\&mu;}}_{k+1,i}<1\end{array}\right.,i=\mathrm{1,2},...,n---\left(22\right)$

μ _k+1,i＝trace(N _k+1)/trace(M _ii,k+！) (23)

N _k+1＝V _k+1-HQ _k+1H ^T-lR _k+1 (24)

M _k+1＝J _k+1HH ^T (25)

$J_{k+1}=\underset{i=1}{\overset{2n}{\mathrm{\&Sigma;}}}{W}_i^c({\mathrm{\&xi;}}_{k+1|k,i}-{\mathrm{\&xi;}}_{0,k+1|k}){({\mathrm{\&xi;}}_{k+1|k,i}-{\mathrm{\&xi;}}_{0,k+1|k})}^T---\left(26\right)$

$V_{k+1}=\left\{\begin{array}{c}{e}_{k+1}{e}_{k+1}^T,k=0\\ \frac{\mathrm{\&rho;}{V}_k+e_{k+1}{e}_{k+1}^T}{1+\mathrm{\&rho;}},k\&GreaterEqual;1\end{array}\right.---\left(27\right)$

4. measure and upgrade

The square root of state one-step prediction covariance matrix calculates:

$S_{k+1|k}=\mathrm{qr}\left\{\right[\mathrm{diag}\left({\mathrm{\&mu;}}_{k+1}\right)\sqrt{W_1^c}({\mathrm{\&xi;}}_{1:2n,k+1|k}-{\mathrm{\&xi;}}_{0,k+1|k}),G_k\sqrt{Q_k}\left]\right\}---\left(28\right)$

Cross-covariance calculates:

$P_{\mathrm{XZ},k+1}=\left(S_{k+1|k}{S}_{k+1|k}^T\right)H^T---\left(29\right)$

The square root of output covariance matrix calculates:

U _k+1＝HS _k+1|k (30)

$S_{Z,k+1}=\mathrm{qr}\left\{\right[U_{k+1},\sqrt{R_{k+1}}\left]\right\}---\left(31\right)$

State gain matrix calculates:

$K_{k+1}=(P_{\mathrm{XZ},k+1}/S_{Z,k+1}^T)/S_{Z,k+1}---\left(32\right)$

State estimation value is calculated:

${\hat{X}}_{k+1|k+1}={\hat{X}}_{k+1|k}+K_{k+1}{e}_{k+1}---\left(33\right)$

The square root of state error covariance matrix calculates:

$S_{k+1}=S_{k+1|k}-S_{k+1|k}{U}_{k+1}^T{S}_{z,k+1}^{-T}{(S_{z,k+1}+\sqrt{R_{k+1}})}^{-1}{U}_{k+1}---\left(34\right)$

The calculation of parameter of using in said process is as follows:

$W_0^m=\frac{\mathrm{\&lambda;}}{n+\mathrm{\&lambda;}}---\left(35\right)$

$W_0^c=W_0^m+(1-{\mathrm{\&alpha;}}^2+\mathrm{\&beta;})---\left(36\right)$

$W_i^c=W_i^m=\frac{1}{2(n+\mathrm{\&lambda;})}---\left(37\right).$

7. a kind of Satellite attitude determination method based on improved self-adaptation square root UKF algorithm according to claim 6, is characterized in that the attitude algorithm method described in step 8 is:

θ＝arcsin(-2(q ₁q ₃-q ₂q ₄)) (39)

$\mathrm{\&psi;}=\arctan (2(q_2{q}_1+q_3{q}_4)/(1-2(q_3^2+q_2^2))---\left(40\right)$

In formula, θ, ψ represents respectively roll angle, the angle of pitch and crab angle.