CN101846510A - High-precision satellite attitude determination method based on star sensor and gyroscope - Google Patents

High-precision satellite attitude determination method based on star sensor and gyroscope Download PDF

Info

Publication number
CN101846510A
CN101846510A CN 201010194288 CN201010194288A CN101846510A CN 101846510 A CN101846510 A CN 101846510A CN 201010194288 CN201010194288 CN 201010194288 CN 201010194288 A CN201010194288 A CN 201010194288A CN 101846510 A CN101846510 A CN 101846510A
Authority
CN
China
Prior art keywords
overbar
omega
satellite
formula
centerdot
Prior art date
Application number
CN 201010194288
Other languages
Chinese (zh)
Other versions
CN101846510B (en
Inventor
杨静
魏明坤
Original Assignee
北京航空航天大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 北京航空航天大学 filed Critical 北京航空航天大学
Priority to CN 201010194288 priority Critical patent/CN101846510B/en
Publication of CN101846510A publication Critical patent/CN101846510A/en
Application granted granted Critical
Publication of CN101846510B publication Critical patent/CN101846510B/en

Links

Abstract

The present invention discloses a high-precision satellite attitude determination method based on star sensor and gyroscope, which comprises the following steps: step 1. establishing a status equation of a satellite attitude determination system; step 2. establishing a measurement equation of the satellite attitude determination system; step 3. performing an online real-time model error estimation through predictive filtering; and step 4. performing a status estimation on a compensated model through 2-order interpolation filtering to obtain the attitude of a satellite. By applying predictive filtering to performing an online real-time model error estimation and correcting the system model, the invention overcomes the shortcoming existing in the conventional estimation process that error is processed into zero-mean white noise; and in addition, the invention can process any nonlinear system and noise conditions to obtain an estimation result of higher precision and is applicable to the field of high-precision attitude determination.

Description

A kind of high-precision satellite attitude based on star sensor and gyro is determined method

Technical field

The present invention relates to a kind of high-precision satellite attitude and determine method, belong to the spacecraft high-precision attitude and determine technical field based on star sensor and gyro.

Background technology

Attitude of flight vehicle determines that system is the important component part in the attitude of flight vehicle control system, and its precision and reliability are the important assurance and the prerequisites of the long-term normal operation of aircraft.This system mainly determines that by attitude sensor and corresponding attitude algorithm forms, and for the influence to attitude determination accuracy of the uncertainty that overcomes reference vector, the method for state estimation determines that at high-precision attitude occasion has obtained using widely.

Method at present commonly used be EKF (Extended Kalman Filter, EKF), EKF has fast convergence rate, algorithm is simple, can real-time resolving etc. advantage.But EKF is because adopted nonlinear model Taylor has been launched approximate method, non-linear when strong estimated bias bigger.And EKF need calculate Jacobian matrix, requires the function must can be little, has limited its usable range.In addition, EKF is treated to the zero-mean white noise with system noise, and this just requires the modeling must be very accurate, if the noise in the real system is not a white Gaussian noise, handles and be used as white noise, then will cause estimated accuracy seriously to descend.Nonlinear prediction filtering (Nonlinear Predict Filter, NPF) system model error or the unknown influence of importing have been considered, has estimated capacity to model error, compare with other non-linear filtering methods such as EKF, predictive filtering can On-line Estimation go out unknown model error, overcome the limitation that model error and system noise is made as white Gaussian noise, in attitude of flight vehicle estimation field, obtained application.But the predictive filtering speed of convergence is relatively slow, and the inaccurate meeting of the model error that estimates during convergence impacts filtering convergence and precision.

Summary of the invention

The objective of the invention is to determine the aspect the deficiencies in the prior art, propose a kind of high-precision satellite attitude and determine method based on star sensor and gyro in order to solve the aircraft high-precision attitude.

A kind of high-precision satellite attitude based on star sensor and gyro of the present invention is determined method, comprises following step:

Step 1: the state equation of setting up Satellite Attitude Determination System;

Step 2: the measurement equation of setting up Satellite Attitude Determination System;

Step 3: utilize predictive filtering online in real time estimation model error;

Step 4: utilize second order interpolation filtering to carry out state estimation to the model after the compensation, obtain the attitude of satellite.

The invention has the advantages that:

(1) adopted predictive filtering online in real time estimation model error and update the system model, having overcome in the conventional estimated process Error processing is the shortcoming of zero-mean white noise.

(2) adopted the second order interpolation filtering algorithm, be a kind of polynomial interpolation filtering algorithm higher than EKF precision, and do not need to calculate partial derivative, the scope of application is wider;

(3) the present invention has adopted predictive filtering to estimate and the bucking-out system model error, can handle any nonlinear system noise situations of unifying, and obtains more high-precision estimated result, is applicable to that high-precision attitude determines the field.

Description of drawings

Fig. 1 is a method flow diagram of the present invention;

Fig. 2 a is the angle of pitch graph of errors comparison diagram of method of the present invention and EKF method simulation result;

Fig. 2 b is the crab angle graph of errors comparison diagram of method of the present invention and EKF method simulation result;

Fig. 2 c is the roll angle graph of errors comparison diagram of method of the present invention and EKF method simulation result;

Fig. 3 a is the angle of pitch graph of errors comparison diagram of method of the present invention and second order interpolation filtering method simulation result;

Fig. 3 b is the crab angle graph of errors comparison diagram of method of the present invention and second order interpolation filtering method simulation result;

Fig. 3 c is the roll angle graph of errors comparison diagram of method of the present invention and second order interpolation filtering method simulation result.

Embodiment

The present invention is described in further detail below in conjunction with drawings and Examples.

The present invention is that a kind of high-precision satellite attitude based on starlight and gyro is determined method, utilize gyro and star sensor as attitude sensor, utilize the thought of predictive filtering to estimate the model error of Satellite Attitude Determination System in real time and compensate in real time, utilize the higher second order interpolation Filtering Estimation of precision to go out the attitude quaternion of satellite body coordinate system again with respect to inertial coordinates system, obtain the attitude of satellite through resolving then, flow process comprises following step as shown in Figure 1:

Step 1: the state equation of setting up Satellite Attitude Determination System;

The attitude of coming together to describe satellite with attitude dynamic equations and kinematical equation changes, and attitude dynamic equations can derive from the momentum moment formula and the theorem of rigid body, i.e. Euler-Newton method: according to the rigid body moment of momentum theorem, the attitude dynamic equations of satellite is:

J ω · + ω × Jω = G gb + G mb + G ab + G sb - - - ( 1 )

In the formula, ω=[ω x, ω y, ω z] TBe the rotational angular velocity of satellite body coordinate system to inertial coordinates system, ω x, ω y, ω zFor the satellite body coordinate system ties up to satellite body coordinate system lower edge x to inertial coordinate, y, the component of three axles of z; J is the inertial tensor matrix of satellite; G MbBe geomagnetic torque; G AbBe aerodynamic moment; G SbBe the solar pressure square; G GbBe gravity gradient torque;

G GbBe specially:

G gb = 3 μ r 5 R b × JR b - - - ( 2 )

In the formula, μ is a geocentric gravitational constant; R is the earth's core distance of satellite; R bBe the position vector of satellite under the satellite body coordinate system, be specially:

R b = C i b R i - - - ( 3 )

In the formula, R iBe the position vector of satellite under inertial coordinates system, For be tied to the transition matrix of satellite body coordinate system from inertial coordinate, be specially:

The earth's core distance True anomaly Eccentric anomaly Mean anomaly M=n (t-t0); N is the flat movement velocity of satellite; t 0For satellite arrives the perigean time; T is the time; Be the transition matrix of geocentric orbital reference system, be specially to inertial coordinates system:

C o i = cos Ω o cos w o - sin w o cos i o sin Ω o - cos Ω o sin w o - cos w o cos i o sin Ω o sin Ω o sin i o sin Ω o cos w o + sin w o cos i o cos Ω o - sin Ω o sin w o + cos w o cos i o cos Ω o - cos Ω o sin i o sin w o sin i o sin w o cos i o cos i o

Wherein, Ω oBe the ascending node of satellite orbit right ascension; i oBe inclination of satellite orbit; w oBe argument of perigee of satellite orbit.

Because G GbThan geomagnetic torque G Mb, aerodynamic moment G AbWith solar pressure square G SbInfluence to Satellite Attitude Determination System is much bigger, so with G Mb, G AbAnd G SbWith resultant moment d GExpression, as the model error of modeling not, then the attitude of satellite kinetics equation of formula (1) is written as:

ω · = J - 1 ( - ω × Jω + G gb ) + J - 1 d G - - - ( 4 )

It is the equation that attitude parameter changes in the attitude maneuver process that attitude motion is learned equation, because it is linear differential equation that the attitude motion that hypercomplex number is represented is learned equation, do not contain trigonometric function, no singular point problem, so the Satellite Attitude Movement equation that adopts hypercomplex number to describe is as follows:

q · = 1 2 Ω ( ω ) q = 1 2 Ξ ( q ) ω - - - ( 5 )

In the formula, Be the attitude quaternion of satellite body coordinate system with respect to inertial coordinates system; q 13=[q 1q 2q 3] T, I 3 * 3Be unit matrix; [ω *], [q 13*] represent respectively by vectorial ω, q 13The antisymmetric matrix that constitutes of component,

The state equation that is made of Satellite Attitude Determination System formula (4) and formula (5) is:

x · = J - 1 ( - ω × Jω + G gb ) 1 2 Ξ ( q ) ω + G · d G + W - - - ( 6 )

In the formula: x is a state vector, x=[ω T, q T] TG is the error perturbation matrix, System noise w (k)~(0, Q (k)), promptly w (k) obeys that average is 0, variance is the Gaussian distribution of Q (k).

With formula (6) brief note be:

x · = f ( x , d G , w ) - - - ( 7 )

Step 2: the measurement equation of setting up Satellite Attitude Determination System;

Adopt three gyros and star sensor as attitude sensor among the present invention.

1) gyro;

The measurement equation of gyro is:

g 1 ( x ( k ) ) = ω x ω y ω z + v x v y v z = g ω ( x ( k ) ) + v ω ( k ) - - - ( 8 )

In the formula, g 1Be the measurement output of gyro, note (i=x, y z) are the rotational angular velocity of satellite body coordinate system to inertial coordinates system, v to ω i i(i=x, y z) are the white Gaussian noise of gyro to measure.

(2) star sensor;

The present invention adopts two star sensors to measure, and the unit direction vector of the primary optical axis of two star sensitivities is respectively l in inertial coordinates system I1And l I2, then the measurement equation of star sensor is:

g 2 = C i b l i 1 + v s 1 = g s 1 ( x ( k ) ) + v s 1 g 3 = C i b l i 2 + v s 2 = g s 2 ( x ( k ) ) + v s 2 - - - ( 9 )

In the formula, v S1, v S2Be the noise that star sensor is measured, mainly be thought of as white Gaussian noise here; g 1, g 2Be respectively the measurement output of two star sensors; Note Be the transition matrix that is tied to the satellite body coordinate system by inertial coordinate, available attitude quaternion is represented, is specially:

C i b = q 1 2 - q 2 2 - q 3 2 + q 4 2 2 ( q 1 q 2 + q 3 q 4 ) 2 ( q 1 q 3 - q 2 q 4 ) 2 ( q 1 q 2 - q 3 q 4 ) - q 1 2 + q 2 2 - q 3 2 + q 4 2 2 ( q 2 q 3 + q 1 q 4 ) 2 ( q 1 q 3 + q 2 q 4 ) 2 ( q 2 q 3 - q 1 q 4 ) - q 1 2 - q 2 2 + q 3 2 + q 4 2 - - - ( 10 )

(3) observation equation of setting up Satellite Attitude Determination System is:

y ( k ) = g ω ( x ( k ) ) g s 1 ( x ( k ) ) g s 2 ( x ( k ) ) + v ω ( k ) v s 1 ( k ) v s 2 ( k ) - - - ( 11 )

With formula (11) brief note be:

y(k)=g(x(k),v(k)) (12)

Wherein, Be that v (k) obedience average is 0, variance is the Gaussian distribution of R (k), and w (k), and v (k) is separate.

Step 3: utilize predictive filtering online in real time estimation model error;

Model error Estimator be:

d ^ G ( k ) = - ( { Λ ( T ) U [ x ^ ( k ) ] } T R - 1 { Λ ( T ) U [ x ^ ( k ) ] } + W ) - 1 { Λ ( T ) U [ x ^ ( k ) ] } T R - 1 (13)

x { y ^ ( k ) + z [ x ^ ( k ) , T ] - y ( k + 1 ) }

In the formula: Estimate that for the output of gyro and two star sensors T is the filtering cycle, U is a sensitivity matrix, and Λ (T) is a diagonal matrix, It is a column vector.W is a weighting matrix, preestablishes numerical value.

Wherein, Be specially:

y ^ ( k ) = g [ x ^ ( k ) v ‾ ( k ) ] - - - ( 14 )

In the formula: Be the estimated value of k quantity of state constantly, The average of expression v (k).

Sensitivity matrix U is:

U = L G [ L f 0 ( g w ) ] L G [ L f 1 ( g s 1 ) ] L G [ L f 1 ( g s 2 ) ] , - - - ( 15 )

Being specially of relevant Lie derivative:

L G [ L f 0 ( g w ) ] = J - 1

L G [ L 1 f ( g s 1 ) ] = ∂ L 1 f ( g s 1 ) ∂ x ^ · G = - Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) J - 1

L G [ L 1 f ( g s 2 ) ] = ∂ L 1 f ( g s 2 ) ∂ x ^ · G = - Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) J - 1

Diagonal matrix Λ (T) is:

Λ ( T ) = Λ w 0 0 0 Λ s 1 0 0 0 Λ s 2 - - - ( 16 )

Wherein: Λ w=TI 3 * 3, Λ s 1 = T 2 2 · I 3 × 3 , Λ s 2 = T 2 2 · I 3 × 3 , ;

Column vector For:

z [ x ^ ( k ) , T ] = z ω z s 1 z s 2 = T · L 1 f ( g ω ) T · L 1 f ( g s 1 ) + T 2 2 L 2 f ( g s 1 ) T · L 1 f ( g s 2 ) + T 2 2 · L 2 f ( g s 2 ) - - - ( 17 )

Each rank Lie derivative is in the following formula:

L 1 f ( g ω ) = J - 1 [ - ω ^ × J ω ^ + G gb ]

L 1 f ( g s 1 ) = - Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) ω ^

L 2 f ( g s 1 ) = [ ω ^ × ] · Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) ω ^ - Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) × J - 1 [ - ω ^ × J ω ^ + G gb ]

L 1 f ( g s 2 ) = - Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) ω ^

L 2 f ( g s 2 ) = [ ω ^ × ] · Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) ω ^ - Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) × J - 1 [ - ω ^ × J ω ^ + G gb ]

Wherein: Γ ( a ) = - [ a ] × - a a T 0 .

At last, obtain the estimated value of model error according to formula (13)

Step 4: utilize second order interpolation filtering to carry out state estimation to the model after the compensation, obtain the attitude of satellite;

Filtering interpolation is to utilize the Stirling interpolation formula that nonlinear function is similar to, and is a kind of method of estimation based on minimum mean square error criterion.

Be specially:

With the model error that obtains in the step 3 Substitution state equation (7) compensates, and is write the state equation and the measurement equation of the Satellite Attitude Determination System model after the compensation as discrete non-linear form, for:

x k+1=f(x k,d G,w k)

(18)

y k=g(x k,v k)

If filter value The error variance battle array Square root be Promptly Be Cholesky decompose, the square root of model error variance battle array Q is S w, the square root of measuring the error variance battle array R of noise is S v, the variance battle array of status predication error Square root be That is:

Q = S w S w T , R = S v S v T , (19)

P ‾ = S ‾ x S - T x , P ^ = S ^ x S ^ x T .

The one-step prediction of state and error variance battle array thereof:

x ‾ k + 1 = h 2 - n x - n w h 2 f ( x ^ k , d ^ G , w ‾ k )

+ 1 2 h 2 Σ p = 1 n x ( f ( x ^ k + h s ^ x , p , d ^ G , w ‾ k ) + f ( x ^ k - h s ^ x , p , d ^ G , w ‾ k ) ) - - - ( 20 )

+ 1 2 h 2 Σ p = 1 n v ( f ( x ^ k , d ^ G , w ‾ k + h s ^ w , p ) + f ( x ^ k , d ^ G , w ‾ k - h s ^ w , p ) )

Wherein, The average of expression w (k), n xThe dimension of expression state vector, n wThe dimension of expression system noise will guarantee that estimated result is not have inclined to one side estimation interpolation step-length to get h 2=3.

The Cholesky of status predication error variance battle array decomposes Matrix after the Householder conversion:

S ‾ ( k + 1 ) = [ S x x ^ ( 1 ) ( k ) S xw ( 1 ) ( k ) S x x ^ ( 2 ) ( k ) S xw ( 2 ) ( k ) ] - - - ( 21 )

Wherein,

S x x ^ ( 1 ) ( k ) = { S x x ^ ( 1 ) ( k ) ( i , j ) } = { ( f i ( x ^ k + h s ^ x , j , d ^ G , w ‾ k ) - f i ( x ^ k - h s ^ x , j , d ^ G , w ‾ k ) ) / 2 h } ,

S x w ( 1 ) ( k ) = { S x w ( 1 ) ( k ) ( i , j ) } = { ( f i ( x ^ k , d ^ G , w ‾ k + h s ^ w , j ) - f i ( x ^ k , d ^ G , w ‾ k - h s ^ w , j ) ) / 2 h } ,

S x x ^ ( 2 ) ( k ) = { S x x ^ ( 2 ) ( k ) ( i , j ) } = { h 2 - 1 2 h 2 ( f i ( x ^ k + h s ^ x , j , d ^ G , w ‾ k )

+ f i ( x ^ k - h s ^ x , j , d ^ G , w ‾ k ) - 2 f i ( x ^ k , d ^ G , w ‾ k ) ) } ,

S x w ( 2 ) ( k ) = { S x w ( 2 ) ( k ) ( i , j ) } = { h 2 - 1 2 h 2 ( f i ( x ^ k , d ^ G , w ‾ k + h s ^ w , j )

+ f i ( x ^ k , d ^ G , w ‾ k - h s ^ w , j ) - 2 f i ( x ^ k , d ^ G , w ‾ k ) ) } ,

Wherein, Expression J row, Expression J row, Expression S wJ row.The one-step prediction of observed reading is:

y ‾ k + 1 = h 2 - n x - n v h 2 g ( x ‾ k + 1 , v ‾ k + 1 )

+ 1 2 h 2 Σ p = 1 n x ( g ( x ‾ k + 1 + h s ‾ x , p , v ‾ k + 1 ) + g ( x ‾ k + 1 - h s ‾ x , p , v ‾ k + 1 ) ) - - - ( 22 )

+ 1 2 h 2 Σ p = 1 n v ( g ( x ‾ k + 1 , v ‾ k + 1 + h s ‾ v , p ) + g ( x ‾ k + 1 , v ‾ k + 1 - h s ‾ v , p ) )

Wherein, n vBe to measure the noise dimension. The Cholesky of error variance battle array to decompose be S y(k+1) matrix after the Householder conversion:

S y ( k + 1 ) = [ S y x ‾ ( 1 ) ( k + 1 ) S yv ( 1 ) ( k + 1 ) S y x ‾ ( 2 ) ( k + 1 ) S yv ( 2 ) ( k + 1 ) ] - - - ( 23 )

Wherein,

S y x ‾ ( 1 ) ( k + 1 ) = { S y x ‾ ( 1 ) ( k + 1 ) ( i , j ) } = { ( g i ( x ‾ k + 1 + h s ‾ x , j , v ‾ k + 1 ) - g i ( x ‾ k + 1 - h s ‾ x , j , v ‾ k + 1 ) ) / 2 h } ,

S yv ( 1 ) ( k + 1 ) = { S yv ( 1 ) ( k + 1 ) ( i , j ) } = { ( g i ( x ‾ k + 1 , v ‾ k + 1 + h s ‾ v , j ) - g i ( x ‾ k + 1 , v ‾ k + 1 - h s ‾ v , j ) ) / 2 h } ,

S y x ‾ ( 2 ) ( k + 1 ) = { S y x ‾ ( 2 ) ( k + 1 ) ( i , j ) } = { h 2 - 1 2 h 2 ( g i ( x ‾ k + 1 + h s ‾ x , j , v ‾ k + 1 )

+ g i ( x ‾ k + 1 - h s ‾ x , j , v ‾ k + 1 ) - 2 g i ( x ‾ k + 1 , v ‾ k + 1 ) ) } ,

S yv ( 2 ) ( k + 1 ) = { S yv ( 2 ) ( k + 1 ) ( i , j ) } = { h 2 - 1 2 h 2 ( g i ( x ‾ k + 1 , v ‾ k + 1 + h s ‾ v , j )

+ g i ( x ‾ k + 1 , v ‾ k + 1 - h s - v , j ) - 2 g i ( x ‾ k + 1 , v ‾ k + 1 ) ) } ,

In the formula: Expression S vJ row.

The cross covariance battle array of state and measurement is:

P xy ( k + 1 ) = S ‾ x ( k + 1 ) S y x ‾ ( k + 1 ) T - - - ( 24 )

Obtain gain matrix K by formula (24) K+1For:

K k+1=P xy(k+1)[S y(k+1)S y(k+1) T] -1 (25)

Be updated to based on state and the estimation error variance battle array measured:

x ^ k + 1 = x - k + 1 + K k + 1 ( y k + 1 - y ‾ k + 1 ) - - - ( 26 )

P ^ ( k + 1 ) = [ S ‾ x ( k + 1 ) - K k + 1 S y x ‾ ( 1 ) ( k + 1 ) ] [ S ‾ x ( k + 1 ) - K k + 1 S yx ( 1 ) ( k + 1 ) ] T + K k + 1 S yv ( 1 ) ( k + 1 ) [ K k + 1 S yv ( 1 ) ( k + 1 ) ] T

+ K k + 1 S yx ( 2 ) ( k + 1 ) [ K k + 1 S yx ( 2 ) ( k + 1 ) ] T + K k + 1 S yv ( 2 ) ( k + 1 ) [ K k + 1 S yv ( 2 ) ( k + 1 ) ] T - - - ( 27 )

Obtain the state estimation value thus Wherein Be the attitude quaternion of satellite, each attitude quaternion carried out normalized after estimating, and resolve according to formula (28) (29) (30) and to obtain attitude angle;

The angle of pitch: θ=-arcsin (2 Q1q3-2 Q2q4) (28)

Crab angle: Ψ = arctan ( 2 q 2 q 3 - 2 q 1 q 4 q 3 2 + q 4 2 - q 1 2 - q 2 2 ) - - - ( 29 )

Roll angle:

This method is because there is pair model error ability of estimation in real time, and therefore the compensation model error can obtain more high-precision Satellite Attitude Estimation in real time.

Simultaneously, this method can also be applicable to other similar systems that has model error.

When having model error, under same simulated conditions, the graph of errors of using attitude that the present invention determines satellite and application extension Kalman filtering (EKF) simulation result is shown in Fig. 2 a, Fig. 2 b, Fig. 2 c, and application the present invention determines that the graph of errors of the attitude of satellite and second order interpolation filtering (DD2) simulation result is shown in Fig. 3 a, Fig. 3 b, Fig. 3 c.

Horizontal ordinate all is simulation times among Fig. 2 a and Fig. 3 a, ordinate is the evaluated error curve (actual value and estimated value poor) of the angle of pitch, as can be seen from the figure, the angle of pitch evaluated error of EKF and DD2 is in 10 rads, and evaluated error of the present invention is in 1 rad.

Horizontal ordinate all is simulation times among Fig. 2 b and Fig. 3 b, and ordinate is the evaluated error curve of crab angle; The horizontal ordinate of Fig. 2 c and Fig. 3 c all is simulation times, and ordinate is the evaluated error curve of roll angle; As can be seen from the figure, the crab angle of EKF and DD2 and the concussion of roll angle evaluated error are more violent, even sometimes error surpasses 15 rads, and evaluated error of the present invention is because the estimation model error also compensates in real time, so the evaluated error angle is still in 2 rads; Therefore the method precision of the present invention's proposition is significantly improved.

Claims (1)

1. the high-precision satellite attitude based on starlight and gyro is determined method, it is characterized in that, adopt gyro and star sensor as attitude sensor, the model error that utilizes the predictive filtering algorithm to estimate Satellite Attitude Determination System in real time also compensates in real time, utilize second order interpolation filtering to carry out state estimation again, specifically comprise following step:
Step 1: the state equation of setting up Satellite Attitude Determination System;
The attitude dynamic equations of satellite is:
J ω · + ω × Jω = G gb + G gb + G gb + G sb - - - ( 1 )
In the formula, ω=[ω x, ω y, ω z] TBe the rotational angular velocity of satellite body coordinate system to inertial coordinates system, ω x, ω y, ω zFor the satellite body coordinate system ties up to satellite body coordinate system lower edge x to inertial coordinate, y, the component of three axles of z; J is the inertial tensor matrix of satellite; G MbBe geomagnetic torque; G AbBe aerodynamic moment; G SbBe the solar pressure square; G GbBe gravity gradient torque;
G GbBe specially:
G gb = 3 μ r 5 R b × J R b - - - ( 2 )
In the formula, μ is a geocentric gravitational constant; R is the earth's core distance of satellite; R bBe the position vector of satellite under the satellite body coordinate system, be specially:
R b = C i b R i - - - ( 3 )
In the formula, R iBe the position vector of satellite under inertial coordinates system, For be tied to the transition matrix of satellite body coordinate system from inertial coordinate, be specially:
The earth's core distance True anomaly Eccentric anomaly Mean anomaly M=n (t-t 0); N is the flat movement velocity of satellite; t 0For satellite arrives the perigean time; T is the time; Be the transition matrix of geocentric orbital reference system, be specially to inertial coordinates system:
C o i = cos Ω o cos w o - sin w o cos i o sin Ω o - cos Ω o sin w o - cos w o cos i o sin Ω o sin Ω o sin i o sin Ω o cos w o + sin w o cos i o cos Ω o - sin Ω o sin w o + cos w o cos i o cos Ω o - cos Ω o sin i o sin w o sin i o sin w o cos i o cos i o
Wherein, Ω oBe the ascending node of satellite orbit right ascension; i oBe inclination of satellite orbit; w oBe argument of perigee of satellite orbit;
With G Mb, G AbAnd G SbResultant moment be expressed as d G, as the model error of modeling not, then the attitude of satellite kinetics equation of formula (1) is written as:
ω · = J - 1 ( - ω × Jω + G gb ) + J - 1 d G - - - ( 4 )
The Satellite Attitude Movement equation that adopts hypercomplex number to describe is as follows:
q · = 1 2 Ω ( ω ) q = 1 2 Ξ ( q ) ω - - - ( 5 )
In the formula, Be the attitude quaternion of satellite body coordinate system with respect to inertial coordinates system; q 13=[q 1q 2q 3] T, I 3 * 3Be unit matrix; [ω *], [q 13*] represent respectively by vectorial ω, q 13The antisymmetric matrix that constitutes of component,
The state equation that is made of Satellite Attitude Determination System formula (4) and formula (5) is:
x · = J - 1 ( - ω × Jω + G gb ) 1 2 Ξ ( q ) ω + G · d G + W - - - ( 6 )
In the formula: x is a state vector, x=[ω T, q T] TG is the error perturbation matrix, System noise w (k)~(0, Q (k)), promptly w (k) obeys that average is 0, variance is the Gaussian distribution of Q (k);
With formula (6) brief note be:
x · = f ( x , d G , w ) - - - ( 7 )
Step 2: the measurement equation of setting up Satellite Attitude Determination System;
Adopt three gyros and two star sensors as attitude sensor;
1) gyro;
The measurement equation of gyro is:
g 1 ( x ( k ) ) = ω x ω y ω z + v x v y v z = g ω ( x ( k ) ) + v ω ( k ) - - - ( 8 )
In the formula, g 1Be the measurement output of gyro, note ω i is the rotational angular velocity of satellite body coordinate system to inertial coordinates system, v iBe the white Gaussian noise of gyro to measure, i=x, y, z;
2) star sensor;
The unit direction vector of the primary optical axis of two star sensitivities is respectively l in inertial coordinates system I1And l I2, then the measurement equation of star sensor is:
g 2 = C i b I i 1 + v s 1 = g s 1 ( x ( k ) ) + v s 1 g 3 = C i b I i 2 + v s 2 = g s 2 ( x ( k ) ) + v s 2 - - - ( 9 )
In the formula, v S1, v S2It is the white Gaussian noise that star sensor is measured; g 1, g 2Be respectively the measurement output of two star sensors; Note Be the transition matrix that is tied to the satellite body coordinate system by inertial coordinate, represent, be specially with attitude quaternion:
C i b = q 1 2 - q 2 2 - q 3 2 + q 4 2 2 ( q 1 q 2 + q 3 q 4 ) 2 ( q 1 q 3 - q 2 q 4 ) 2 ( q 1 q 2 - q 3 q 4 ) - q 1 2 + q 2 2 - q 3 2 + q 4 2 2 ( q 2 q 3 + q 1 q 4 ) 2 ( q 1 q 3 + q 2 q 4 ) 2 ( q 2 q 3 - q 1 q 4 ) - q 1 2 - q 2 2 + q 3 2 + q 4 2 - - - ( 10 )
3) observation equation of setting up Satellite Attitude Determination System is:
y ( k ) = g ω ( x ( k ) ) g s 1 ( x ( k ) ) g s 2 ( x ( k ) ) + v ω ( k ) v s 1 ( k ) v s 2 ( k ) - - - ( 11 )
With formula (11) brief note be:
y(k)=g(x(k),v(k)) (12)
Wherein, Be that v (k) obeys that average is 0, variance is the Gaussian distribution of R (k), and w (k), v (k) is separate;
Step 3: utilize predictive filtering online in real time estimation model error;
Model error Estimator be:
d ^ G ( k ) = - ( { Λ ( T ) U [ x ^ ( k ) ] } T R - 1 { Λ ( T ) U [ x ^ ( k ) ] } + W ) - 1 { Λ ( T ) U [ x ^ ( k ) ] } T R - 1 (13)
× { y ^ ( k ) + z [ x ^ ( k ) , T ] - y ( k + 1 ) }
In the formula: Estimate that for the output of gyro and two star sensors T is the filtering cycle, U is a sensitivity matrix, and Λ (T) is a diagonal matrix, It is a column vector; W is a weighting matrix, preestablishes numerical value;
Wherein, Be specially:
y ^ ( k ) = g [ x ^ ( k ) , v ‾ ( k ) ] - - - ( 14 )
In the formula: Estimated value for k quantity of state constantly; The average of expression v (k);
Sensitivity matrix U is:
U = L G [ L 0 f ( g w ) ] L G [ L 1 f ( g s 1 ) ] L G [ L 1 f ( g s 2 ) ] , - - - ( 15 )
Being specially of relevant Lie derivative:
L G [ L 0 f ( g w ) ] = J - 1
L G [ L 1 f ( g s 1 ) ] = ∂ L 1 f ( g s 1 ) ∂ x ^ · G = - Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) J - 1
L G [ L 1 f ( g s 2 ) ] = ∂ L 1 f ( g s 2 ) ∂ x ^ · G = - Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) J - 1
Diagonal matrix Λ (T) is:
Λ ( T ) = Λ w 0 0 0 Λ s 1 0 0 0 Λ s 2 - - - ( 16 )
Wherein: Λ w=TI 3 * 3, Λ s 1 = T 2 2 · I 3 × 3 , Λ s 2 = T 2 2 · I 3 × 3 , ; Column vector For:
z [ x ^ ( k ) , T ] = z ω z s 1 z s 2 = T · L 1 f ( g ω ) T · L 1 f ( g s 1 ) + T 2 2 · L 2 f ( g s 1 ) T · L 1 f ( g s 2 ) + T 2 2 · L 2 f ( g s 2 ) - - - ( 17 )
Each rank Lie derivative is in the following formula:
L 1 f ( g ω ) = J - 1 [ - ω ^ × J ω ^ + G gb ]
L 1 f ( g s 1 ) = - Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) ω ^
L 2 f ( g s 1 ) = [ ω ^ × ] · Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) ω ^ - Ξ T ( q ^ ) Γ ( l i 1 ) Ξ ( q ^ ) × J - 1 [ - ω ^ × J ω ^ + G gb ]
L 1 f ( g s 2 ) = - Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) ω ^
L 2 f ( g s 2 ) = [ ω ^ × ] · Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) ω ^ - Ξ T ( q ^ ) Γ ( l i 2 ) Ξ ( q ^ ) × J - 1 [ - ω ^ × J ω ^ + G gb ]
Wherein: Γ ( a ) = - [ a × ] - a a T 0 ;
At last, obtain the estimated value of model error according to formula (13)
Step 4: utilize second order interpolation filtering to carry out state estimation to the model after the compensation, obtain the attitude of satellite;
Be specially:
With the model error that obtains in the step 3 Substitution state equation (7) compensates, and is write the state equation and the measurement equation of the Satellite Attitude Determination System model after the compensation as discrete non-linear form, for:
x k+1=f(x k,d G,w k) (18)
y k=g(x k,v k)
If filter value The error variance battle array Square root be Promptly Be Cholesky decompose, the square root of model error variance battle array Q is S w, the square root of measuring the variance battle array R of noise is S v, the variance battle array of status predication error Square root be That is:
Q = S w S w T , R = S v S v T , (19)
P ‾ = S ‾ x S ‾ x T , P ^ = S ^ x S ^ x T .
The one-step prediction of state and error variance battle array thereof:
x ‾ k + 1 = h 2 - n x - n w h 2 f ( x ^ k , d ^ G , w ‾ k )
+ 1 2 h 2 Σ p = 1 n x ( f ( x ^ k + h s ^ x , p , d ^ G , w ‾ k ) + f ( x ^ k - h s ^ x , p , d ^ G , w ‾ k ) ) - - - ( 20 )
+ 1 2 h 2 Σ p = 1 n v ( f ( x ^ k , d ^ G , w ‾ k + h s ^ w , p ) + f ( x ^ k , d ^ G , w ‾ k - h s ^ w , p ) )
Wherein, The average of expression w (k), n xThe dimension of expression state vector, n wThe dimension of expression system noise, interpolation step-length h 2=3;
The Cholesky of status predication error variance battle array decomposes Matrix after the Householder conversion:
S _ ( k + 1 ) = [ S x x ^ ( 1 ) ( k ) S xw ( 1 ) ( k ) S x x ^ ( 2 ) ( k ) S xw ( 2 ) ( k ) ] - - - ( 21 )
Wherein,
S x x ^ ( 1 ) ( k ) = { S x x ^ ( 1 ) ( k ) ( i , j ) } = { ( f i ( x ^ k + h s ^ x , j , d ^ G , w ‾ k ) - f i ( x ^ k - h s ^ x , j , d ^ G , w ‾ k ) ) / 2 h } ,
S xw ( 1 ) ( k ) = { S xw ( 1 ) ( k ) ( i , j ) } = { ( f i ( x ^ k , d ^ G , w ‾ k + h s ^ w , j ) - f i ( x ^ k , d ^ G , w ‾ k - h s ^ w , j ) ) / 2 h } ,
S x x ^ ( 2 ) ( k ) = { S x x ^ ( 2 ) ( k ) ( i , j ) } = { h 2 - 1 2 h 2 ( f i ( x ^ k + h s ^ x , j , d ^ G , w ‾ k )
+ f i ( x ^ k - h s ^ x , j , d ^ G , w ‾ k ) - 2 f i ( x ^ k , d ^ G , w ‾ k ) ) } ,
S x w ( 2 ) ( k ) = { S x w ( 2 ) ( k ) ( i , j ) } = { h 2 - 1 2 h 2 ( f i ( x ^ k , d ^ G , w ‾ k + h s ^ w , j )
+ f i ( x ^ k , d ^ G , w ‾ k - h s ^ w , j ) - 2 f i ( x ^ k , d ^ G , w ‾ k ) ) } ,
Wherein, Expression J row, Expression J row, Expression S wJ row; The one-step prediction of observed reading is:
y ‾ k + 1 = h 2 - n x - n v h 2 g ( x ‾ k + 1 , v ‾ k + 1 )
+ 1 2 h 2 Σ p = 1 n x ( g ( x ‾ k + 1 + h s ‾ x , p , v ‾ k + 1 ) + g ( x ‾ k + 1 - h s ‾ x , p , v ‾ k + 1 ) ) - - - ( 22 )
+ 1 2 h 2 Σ p = 1 n v ( g ( x ‾ k + 1 , v ‾ k + 1 + h s ‾ v , p ) + g ( x ‾ k + 1 , v ‾ k + 1 - h s ‾ v , p ) )
Wherein, n vBe to measure the noise dimension; The Cholesky of error variance battle array to decompose be S y(k+1) matrix after the Householder conversion:
S y ( k + 1 ) = [ S y x ‾ ( 1 ) ( k + 1 ) S yv ( 1 ) ( k + 1 ) S y x ‾ ( 2 ) ( k + 1 ) S yv ( 2 ) ( k + 1 ) ] - - - ( 23 )
Wherein,
S y x ‾ ( 1 ) ( k + 1 ) = { S y x ‾ ( 1 ) ( k + 1 ) ( i , j ) } = { ( g i ( x ‾ k + 1 + h s ‾ x , j , v ‾ k + 1 ) - g i ( x ‾ k + 1 - h s ‾ x , j , v ‾ k + 1 ) ) / 2 h } ,
S yv ( 1 ) ( k + 1 ) = { S yv ( 1 ) ( k + 1 ) ( i , j ) } = { ( g i ( x ‾ k + 1 , v ‾ k + 1 + h s ‾ v , j ) - g i ( x ‾ k + 1 , v ‾ k + 1 - h s ‾ v , j ) ) / 2 h } ,
S y x ‾ ( 2 ) ( k + 1 ) = { S y x ‾ ( 2 ) ( k + 1 ) ( i , j ) } = { h 2 - 1 2 h 2 ( g i ( x ‾ k + 1 + h s ‾ x , j , v ‾ k + 1 )
+ g i ( x ‾ k + 1 - h s ‾ x , j , v ‾ k + 1 ) - 2 g i ( x ‾ k + 1 , v ‾ k + 1 ) ) } ,
S yv ( 2 ) ( k + 1 ) = { S yv ( 2 ) ( k + 1 ) ( i , j ) } = { h 2 - 1 2 h 2 ( g i ( x ‾ k + 1 , v ‾ k + 1 + h s ‾ v , j )
+ g i ( x ‾ k + 1 , v ‾ k + 1 - h s ‾ v , j ) - 2 g i ( x ‾ k + 1 , v ‾ k + 1 ) ) } ,
In the formula: Expression S vJ row;
The cross covariance battle array of state and measurement is:
P xy ( k + 1 ) = S ‾ x ( k + 1 ) S y x ‾ ( k + 1 ) T - - - ( 24 )
Obtain gain matrix K by formula (24) K+1For:
K k+1=P xy(k+1)[S y(k+1)S y(k+1) T] -1 (25)
Be updated to based on state and the estimation error variance battle array measured:
x ^ k + 1 = x - k + 1 + K k + 1 ( y k + 1 - y ‾ k + 1 ) - - - ( 26 )
P ^ ( k + 1 ) = [ S ‾ x ( k + 1 ) - K k + 1 S y x ‾ ( 1 ) ( k + 1 ) ] [ S ‾ x ( k + 1 ) - K k + 1 S yx ( 1 ) ( k + 1 ) ] T + K k + 1 S yv ( 1 ) ( k + 1 ) [ K k + 1 S yv ( 1 ) ( k + 1 ) ] T
+ K k + 1 S yx ( 2 ) ( k + 1 ) [ K k + 1 S yx ( 2 ) ( k + 1 ) ] T + K k + 1 S yv ( 2 ) ( k + 1 ) [ K k + 1 S yv ( 2 ) ( k + 1 ) ] T - - - ( 27 )
Obtain the state estimation value thus Wherein Be the attitude quaternion of satellite, each attitude quaternion carried out normalized after estimating, and resolve according to formula (28) (29) (30) and to obtain attitude angle;
The angle of pitch: θ=-arcsin (2 Q1q3-2 Q2q4) (28)
Crab angle: Ψ = arctan ( 2 q 2 q 3 - 2 q 1 q 4 q 3 2 + q 4 2 - q 1 2 - q 2 2 ) - - - ( 29 )
Roll angle:
CN 201010194288 2010-05-28 2010-05-28 High-precision satellite attitude determination method based on star sensor and gyroscope CN101846510B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN 201010194288 CN101846510B (en) 2010-05-28 2010-05-28 High-precision satellite attitude determination method based on star sensor and gyroscope

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN 201010194288 CN101846510B (en) 2010-05-28 2010-05-28 High-precision satellite attitude determination method based on star sensor and gyroscope

Publications (2)

Publication Number Publication Date
CN101846510A true CN101846510A (en) 2010-09-29
CN101846510B CN101846510B (en) 2013-03-27

Family

ID=42771205

Family Applications (1)

Application Number Title Priority Date Filing Date
CN 201010194288 CN101846510B (en) 2010-05-28 2010-05-28 High-precision satellite attitude determination method based on star sensor and gyroscope

Country Status (1)

Country Link
CN (1) CN101846510B (en)

Cited By (27)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102176159A (en) * 2011-02-28 2011-09-07 哈尔滨工业大学 Satellite attitude control system failure diagnosis device and method based on state observer and equivalent space
CN102288200A (en) * 2011-07-06 2011-12-21 清华大学 Accuracy measurement system for star sensor
CN102288201A (en) * 2011-07-06 2011-12-21 清华大学 Precision measurement method for star sensor
CN102306160A (en) * 2011-07-20 2012-01-04 航天东方红卫星有限公司 Assistant data processing method for improving image positioning precision of stereoscopic plotting camera
CN102436437A (en) * 2011-11-17 2012-05-02 西北工业大学 Quaternion Fourier approximate output method in extreme flight of aircraft based on angular speed
CN102495830A (en) * 2011-11-17 2012-06-13 西北工业大学 Quaternion Hartley approximate output method based on angular velocities for aircraft during extreme flight
CN102495829A (en) * 2011-11-17 2012-06-13 西北工业大学 Quaternion Walsh approximate output method based on angular velocities for aircraft during extreme flight
CN102620886A (en) * 2012-03-27 2012-08-01 南京航空航天大学 Two-step in-orbit recognition rotary inertia estimation method for combined spacecraft
CN102679976A (en) * 2012-05-20 2012-09-19 西安费斯达自动化工程有限公司 Aircraft quaternion revising model based on accelerated speed
CN102865866A (en) * 2012-10-22 2013-01-09 哈尔滨工业大学 Satellite attitude determination method and attitude determination error analytical method based on two star sensors
CN102867130A (en) * 2012-10-11 2013-01-09 西北工业大学 Aircraft modeling method based on variable measurement number maximum information criterion
CN102901977A (en) * 2012-10-24 2013-01-30 北京航天自动控制研究所 Method for determining initial attitude angle of aircraft
CN103267531A (en) * 2013-04-23 2013-08-28 上海卫星工程研究所 Method for high-precision compensation of fiber-optic gyroscope random error
CN103674023A (en) * 2013-12-26 2014-03-26 中国人民解放军国防科学技术大学 Method for dynamically measuring attitude of star sensor based on top accurate angle relevance
CN103954288A (en) * 2014-05-19 2014-07-30 中国人民解放军国防科学技术大学 Determination method for precision responding relation of satellite attitude determination system
CN104020671A (en) * 2014-05-30 2014-09-03 哈尔滨工程大学 Robustness recursion filtering method for aircraft attitude estimation under the condition of measurement interference
CN104034334A (en) * 2014-06-05 2014-09-10 哈尔滨工程大学 Single-star and double-star attitude determination method of small-field star sensor
CN104077364A (en) * 2014-06-11 2014-10-01 上海微小卫星工程中心 Globe limb observation method
CN104833358A (en) * 2015-05-13 2015-08-12 上海交通大学 Method for determining geometric linear gesture of star sensor based on Rodrigues parameters
CN105866459A (en) * 2016-03-25 2016-08-17 中国人民解放军国防科学技术大学 Method for estimating restraining angular velocity of gyro-free inertial measurement system
CN106027904A (en) * 2016-06-29 2016-10-12 中国人民解放军国防科学技术大学 Video satellite attitude motion image compensation method
CN106767767A (en) * 2016-11-23 2017-05-31 上海航天控制技术研究所 A kind of micro-nano multimode star sensor system and its data fusion method
CN106843261A (en) * 2016-10-25 2017-06-13 南京航空航天大学 A kind of tensor product interpolation modeling of morphing aircraft changeover portion and control method
CN107065917A (en) * 2017-06-06 2017-08-18 上海微小卫星工程中心 Near space attitude motion of spacecraft characteristic descriptive model and its modeling method
CN107942090A (en) * 2017-12-28 2018-04-20 北京航空航天大学 A kind of spacecraft Attitude rate estimator method based on fuzzy star chart extraction Optic flow information
CN108507527A (en) * 2018-03-07 2018-09-07 北京电子工程总体研究所 A kind of gimbaled nozzle solving of attitude method
CN109781102A (en) * 2019-01-14 2019-05-21 上海卫星工程研究所 Attitude measurement method and system based on double super platforms

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2705448A1 (en) * 1993-05-14 1994-11-25 Hughes Aircraft Co System and method for determining the attitude of a spacecraft using stellar or terrestrial sensors
CN101082494A (en) * 2007-06-19 2007-12-05 北京航空航天大学 Self boundary marking method based on forecast filtering and UPF spacecraft shading device
CN101196398A (en) * 2007-05-25 2008-06-11 北京航空航天大学 Spacecraft posture confirming method based on Euler-q algorithm and DD2 filtering

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2705448A1 (en) * 1993-05-14 1994-11-25 Hughes Aircraft Co System and method for determining the attitude of a spacecraft using stellar or terrestrial sensors
CN101196398A (en) * 2007-05-25 2008-06-11 北京航空航天大学 Spacecraft posture confirming method based on Euler-q algorithm and DD2 filtering
CN101082494A (en) * 2007-06-19 2007-12-05 北京航空航天大学 Self boundary marking method based on forecast filtering and UPF spacecraft shading device

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
《系统仿真学报》 20100228 冀红霞 等 基于非线性预测滤波的卫星姿态确定 34-38,42 1 第22卷, 2 *

Cited By (45)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102176159B (en) * 2011-02-28 2013-05-01 哈尔滨工业大学 Satellite attitude control system failure diagnosis device and method based on state observer and equivalent space
CN102176159A (en) * 2011-02-28 2011-09-07 哈尔滨工业大学 Satellite attitude control system failure diagnosis device and method based on state observer and equivalent space
CN102288200A (en) * 2011-07-06 2011-12-21 清华大学 Accuracy measurement system for star sensor
CN102288201A (en) * 2011-07-06 2011-12-21 清华大学 Precision measurement method for star sensor
CN102288201B (en) * 2011-07-06 2012-08-22 清华大学 Precision measurement method for star sensor
CN102288200B (en) * 2011-07-06 2012-07-18 清华大学 Accuracy measurement system for star sensor
CN102306160B (en) * 2011-07-20 2013-01-23 航天东方红卫星有限公司 Assistant data processing method for improving image positioning precision of stereoscopic plotting camera
CN102306160A (en) * 2011-07-20 2012-01-04 航天东方红卫星有限公司 Assistant data processing method for improving image positioning precision of stereoscopic plotting camera
CN102495829A (en) * 2011-11-17 2012-06-13 西北工业大学 Quaternion Walsh approximate output method based on angular velocities for aircraft during extreme flight
CN102436437B (en) * 2011-11-17 2014-07-09 西北工业大学 Quaternion Fourier approximate output method in extreme flight of aircraft based on angular speed
CN102436437A (en) * 2011-11-17 2012-05-02 西北工业大学 Quaternion Fourier approximate output method in extreme flight of aircraft based on angular speed
CN102495830A (en) * 2011-11-17 2012-06-13 西北工业大学 Quaternion Hartley approximate output method based on angular velocities for aircraft during extreme flight
CN102620886A (en) * 2012-03-27 2012-08-01 南京航空航天大学 Two-step in-orbit recognition rotary inertia estimation method for combined spacecraft
CN102679976A (en) * 2012-05-20 2012-09-19 西安费斯达自动化工程有限公司 Aircraft quaternion revising model based on accelerated speed
CN102679976B (en) * 2012-05-20 2014-07-09 西安费斯达自动化工程有限公司 Aircraft quaternion revising model based on accelerated speed
CN102867130A (en) * 2012-10-11 2013-01-09 西北工业大学 Aircraft modeling method based on variable measurement number maximum information criterion
CN102867130B (en) * 2012-10-11 2015-04-01 西北工业大学 Aircraft modeling method based on variable measurement number maximum information criterion
CN102865866A (en) * 2012-10-22 2013-01-09 哈尔滨工业大学 Satellite attitude determination method and attitude determination error analytical method based on two star sensors
CN102865866B (en) * 2012-10-22 2015-01-28 哈尔滨工业大学 Satellite attitude determination method and attitude determination error analytical method based on two star sensors
CN102901977A (en) * 2012-10-24 2013-01-30 北京航天自动控制研究所 Method for determining initial attitude angle of aircraft
CN103267531A (en) * 2013-04-23 2013-08-28 上海卫星工程研究所 Method for high-precision compensation of fiber-optic gyroscope random error
CN103267531B (en) * 2013-04-23 2015-11-18 上海卫星工程研究所 A kind of high-accuracy compensation method of optical fiber gyro random error
CN103674023A (en) * 2013-12-26 2014-03-26 中国人民解放军国防科学技术大学 Method for dynamically measuring attitude of star sensor based on top accurate angle relevance
US9316716B2 (en) 2013-12-26 2016-04-19 National University Of Defense Technology Dynamic attitude measurement method of star sensor based on gyro's precise angular correlation
CN103674023B (en) * 2013-12-26 2014-10-15 中国人民解放军国防科学技术大学 Method for dynamically measuring attitude of star sensor based on top accurate angle relevance
CN103954288A (en) * 2014-05-19 2014-07-30 中国人民解放军国防科学技术大学 Determination method for precision responding relation of satellite attitude determination system
CN103954288B (en) * 2014-05-19 2016-08-17 中国人民解放军国防科学技术大学 A kind of Satellite Attitude Determination System precision response relation determines method
CN104020671B (en) * 2014-05-30 2017-01-11 哈尔滨工程大学 Robustness recursion filtering method for aircraft attitude estimation under the condition of measurement interference
CN104020671A (en) * 2014-05-30 2014-09-03 哈尔滨工程大学 Robustness recursion filtering method for aircraft attitude estimation under the condition of measurement interference
CN104034334B (en) * 2014-06-05 2016-09-14 哈尔滨工程大学 Single star of a kind of small field of view star sensor and double star method for determining posture
CN104034334A (en) * 2014-06-05 2014-09-10 哈尔滨工程大学 Single-star and double-star attitude determination method of small-field star sensor
CN104077364A (en) * 2014-06-11 2014-10-01 上海微小卫星工程中心 Globe limb observation method
CN104077364B (en) * 2014-06-11 2017-04-19 上海微小卫星工程中心 Globe limb observation method
CN104833358A (en) * 2015-05-13 2015-08-12 上海交通大学 Method for determining geometric linear gesture of star sensor based on Rodrigues parameters
CN105866459B (en) * 2016-03-25 2018-10-26 中国人民解放军国防科学技术大学 Non-gyro inertial measurement system restriction Attitude rate estimator method
CN105866459A (en) * 2016-03-25 2016-08-17 中国人民解放军国防科学技术大学 Method for estimating restraining angular velocity of gyro-free inertial measurement system
CN106027904A (en) * 2016-06-29 2016-10-12 中国人民解放军国防科学技术大学 Video satellite attitude motion image compensation method
CN106843261A (en) * 2016-10-25 2017-06-13 南京航空航天大学 A kind of tensor product interpolation modeling of morphing aircraft changeover portion and control method
CN106767767A (en) * 2016-11-23 2017-05-31 上海航天控制技术研究所 A kind of micro-nano multimode star sensor system and its data fusion method
CN107065917A (en) * 2017-06-06 2017-08-18 上海微小卫星工程中心 Near space attitude motion of spacecraft characteristic descriptive model and its modeling method
CN107065917B (en) * 2017-06-06 2020-03-17 上海微小卫星工程中心 Near space spacecraft attitude motion characteristic description model and modeling method thereof
CN107942090A (en) * 2017-12-28 2018-04-20 北京航空航天大学 A kind of spacecraft Attitude rate estimator method based on fuzzy star chart extraction Optic flow information
CN107942090B (en) * 2017-12-28 2019-10-29 北京航空航天大学 A kind of spacecraft Attitude rate estimator method for extracting Optic flow information based on fuzzy star chart
CN108507527A (en) * 2018-03-07 2018-09-07 北京电子工程总体研究所 A kind of gimbaled nozzle solving of attitude method
CN109781102A (en) * 2019-01-14 2019-05-21 上海卫星工程研究所 Attitude measurement method and system based on double super platforms

Also Published As

Publication number Publication date
CN101846510B (en) 2013-03-27

Similar Documents

Publication Publication Date Title
Ahmed et al. Accurate attitude estimation of a moving land vehicle using low-cost MEMS IMU sensors
Li et al. Vision-aided inertial navigation for pinpoint planetary landing
CN106289246B (en) A kind of flexible link arm measure method based on position and orientation measurement system
CN103913181B (en) A kind of airborne distributed POS Transfer Alignments based on parameter identification
US20130138264A1 (en) Automotive navigation system and method to utilize internal geometry of sensor position with respect to rear wheel axis
Helmick et al. Path following using visual odometry for a mars rover in high-slip environments
Jiancheng et al. Study on innovation adaptive EKF for in-flight alignment of airborne POS
CN104215259B (en) A kind of ins error bearing calibration based on earth magnetism modulus gradient and particle filter
CN103245359B (en) A kind of inertial sensor fixed error real-time calibration method in inertial navigation system
CN102508278B (en) Adaptive filtering method based on observation noise covariance matrix estimation
CN101788296B (en) SINS/CNS deep integrated navigation system and realization method thereof
CN102692225B (en) Attitude heading reference system for low-cost small unmanned aerial vehicle
Sukkarieh et al. A high integrity IMU/GPS navigation loop for autonomous land vehicle applications
CN100587641C (en) A kind of attitude determination system that is applicable to the arbitrary motion mini system
CN100516775C (en) Method for determining initial status of strapdown inertial navigation system
CN101726295B (en) Unscented Kalman filter-based method for tracking inertial pose according to acceleration compensation
Bonnifait et al. Data fusion of four ABS sensors and GPS for an enhanced localization of car-like vehicles
Ojeda et al. FLEXnav: Fuzzy logic expert rule-based position estimation for mobile robots on rugged terrain
CN104655152B (en) A kind of real-time Transfer Alignments of airborne distributed POS based on federated filter
CN100559125C (en) A kind of spacecraft attitude based on Euler-q algorithm and DD2 filtering is determined method
Mulder et al. Non-linear aircraft flight path reconstruction review and new advances
Han et al. A novel method to integrate IMU and magnetometers in attitude and heading reference systems
CN103076015B (en) A kind of SINS/CNS integrated navigation system based on optimum correction comprehensively and air navigation aid thereof
CN101915580A (en) Self-adaptation three-dimensional attitude positioning method based on microinertia and geomagnetic technology
CN102116628B (en) High-precision navigation method for landed or attached deep sky celestial body detector

Legal Events

Date Code Title Description
PB01 Publication
C06 Publication
SE01 Entry into force of request for substantive examination
C10 Entry into substantive examination
GR01 Patent grant
C14 Grant of patent or utility model
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20130327

Termination date: 20140528

C17 Cessation of patent right