On the Attitude Estimation of Earth Observation Satellites

-wo(Iy- Iz)/I x

c

-wo (I x - I y ) /1 z

g

-3w.2(I - I )/1 y z x

Ix Iy I z are moments of inertia with re-pect to the pr i ncipal axez, b b 1

are constant disturb a nces, and u

2

are control moments.

According to Prof. Brys on and Kortum (1 971), we adopt two system models of a ttitude estima t ion.

The system output i s horizon sensor (HS) measurement (2)

Model-l. The sy s tem dyna mics are described by the Eu l ar's and kinematic equations of motion in the linearized form

where v is a white noise process. It is easy to see th e system (1)-(2) is not completely observable. If we introduce new state variables as they are (Bryson and Kortum, 1971; Yang and Suing, 1980)

p

Xl'

= [4>,

U

= [ u , u

1

N ' 1 --: ...

1

2

u

1

= [ n1,

where y, =

--

2 15

Yl' Y2' r, 2

bd

]

n 2 , n , n41

3

'f - b.llwo c and y,2 = p

+ b2 /

c,

S. H. Yu and Q. P. Chu

2 16

then we can describe the state equation of X, in the form X,

A,X,+ B,U, + F, N,

(3 )

Zl

Cl X, + v

(4 )

r

+ n

r

r

o

As concerning the system measurement,

where 0

A1 =

0

Wo

0

0

-w" 0

0

0

F,=

0

,

0

I)

0

0

g

0

0

a

1

0

0

Cl

0

c

0

0

0

0

0

1

0

0

0

0

0

0

0

0

0

:

B1 '-

[

C

[1

0

:]

0

0

0

0

0

0

oJ

A1 'C 1 " · ••• (A 1 ,)4 C1 ']=5

rank [Cl"

So, the system of Eqs. (3)-(4) is completely observable in the deterministic sense. The discrete-time equations of the system can be given in the following form (Anderson and Moore, 1979) X'k+' = A1k X1k + B1kU'k + F'k N1k Zlk = C'k Xlk + v k

1k

+ p + n

Wolf'

~

-woq> +

p

1

r + n

2

ar + g4> + u, + b, + n3

r

cp + u p

- ci d

p p

2

+ b + n

p

2

+ n

C2X2 + V2

4

(7)

(9)

('0)

where X

2

'=

[ep,

\ji,

U 2 ' = [ u"

=[

p, r, d p ' d r , bp'

u2]

b"

b2]

[ z"

z2' z3]

N2 '= [ n 1 , n 2 , n , n 4 • n p ' nrJ 3 V 2 '=[v"

v 2 • v3J

o w"

1

-w"

0

0

g

0 0

0 c

o

a 0

o 0 o 0 1 0 o , o 0 F2 o 0 o 0 o 0,

0

0

0

0

0

1

0

0

-d

0

r 0 0

0

0

0

0

0

0

00,1

P

o o o

o oo o

l

0

-d

Model-2. For the same system of motion we can construct another model. The model-2 differs from model-' by including the gyro drift parameters d d r b b into the system dynamics p p r and the gyro outputs into the system measurement as shown below

(8 )

2

Z2

o

cov(v ' v ) = R 10 kl k l

+ v

A2 X2 + B2 U2 + B2 b + F2N2

Z2 '=

G,D kl

1l ) =

P

X2

(6)

Po as it s mean and c ovarianc e re s pec-

+ b

p

The system of Eqs. (7)-(8) is complete l y observable, if the constant d isturbances b, and b are considered as 2 known. In order to show this is true we convert the system into state space variable form

b'

tively. The random vectors Nlk and v k are assumed to be white noises with covariances as follows , N

p + d

(5)

with a n initial state, having Xo and

cov(N

now we have the gyro measurements of angular rates in additi on to the previous HS measurement of

0

now the matrix pair (A 1 ' , C1 ,) can be tested to have full rank, i.e.

d

-cl. d

0

0

r1

0 0

001

0

o

0 0 , 000 0 o 000

o o

0

0 0 0

:0

0 1

,

0 I

0 0 0 000

o o o o

0 0 0 0

1

0

o o o

1

0 0

Attitude Estimation of Earth Observation Satellites "-

Thp rank of matrix pair (A 2 ', C2 ') is equal to 8. That is to say the observability condition is satisfied. The discrete-time equations of model-2 are given as follows X2k + 1 = A2k X2k + B 2k U2k + B2k b k ( 11) + F2kN2k Z2k+l

X2k / k = llse of X 2k given (Z2i)' i = O,l •.•• ,k

Then one-step prediction of the state

" + / , is determivector X2k ' i.e. X 2k 1 k ned by the time-updating e quation "

A

X2k + 1 / k = A2kX2k/k + B2k U2k

( 12)

C2k+1X2k+l + V2k + 1

+ B

with E

[x 2 (O)] =

Er E

N21 ]

b

2k k

The subdivision is an iterative (approximation) procedure and has to deal with Eq. (14) only. Suppose we have done nth approximation, now we want to do (n+1)th approximation by assuming

X20

X (O)-X , X2 (O)-X 20 ] 20 2

l N2k ,

217

Gi)kl

ESTIM ATI ON ALGORITHM An estimation algorithm is proposed here, which consists of two parts. Foregoing part. Here the model-l is used to calculate the estimation of state vector Xl by applying linear kal-

The numbers within parentheses indicate to which model the quantities are

man filter theory in the form (Anderson and Moore, 1979)

(1) and 1k Y1k(1) are obtained in the foregoing

X1k + 1 / k

part at time k; 'P ~+1/k(2) is obtained after nth approximation from the subdivision part at time k. We give the nth and (n+1)th approximations of Eq. (14) in more detail form as follows

= AlkXlk/k + B1k U1k

Xl k+l /k+l;;: Xl k+l /k + Kl k+l (Zl k+l - C1k+1Xlk+l/k) K1k + 1

Plk+l/kClk+l(Clk+1Plk+l/kCik+1 + R

+ )-1 1k 1

referred. For example, b"

"n X2k + 1 / k

" " A2kX2k/k + B2k U2k + B2kfk "n 1 ( 1 5) + HX2~+1/k

"n+1 X2k + 1 / k

" " A2kX2k/k + B2k U2k + B2kfk "n ( 1 6) + HX 2k + 1 / k

(13)

P 1k + 1 / k = AlkPlk/kAik + F 1k G1 F ik Plk+l/k+l=Plk+1/k-Klk+1Clk+1P1k+l/k

'"

X10 /- 1 = XO;

P 10 /- 1 = Po

where

In this part we can obtain the esti1\ mates Y1 and "b , which are needed in 1 subdivision part of the algorithm.

Subdivision part.

The main purpose of

this part is to subdivide ~ from their combination

Y,

and ~2

= ~ -

b2 /Wo c

which is found in the foregoing part, and at the same time to obtain the estimates of all of the system variables. Assume we have found a linear least square estimate (llse) of X given 2k Z20' Z21' Z22'···' Z2k' Le.

fk

,

1

b k( 1) 1

= [ -Woe" k(' ) 0

0

0

0

0

0

WoC

0

0

0

0

0

H = B2k [ :

:]

and just for simplicity we take 1\

0

X2k + 1 / k

1\

= X2k / k

Subtracting Eq. (15) from Eq. (16) we obtain "n+1 "n "n "n-1 X2k+1/k-X2k+1/k=H(X2k+1/k-X2k+1/k) ( 17)

S. H. Yu and Q. P. Chu

Repeating this relation m times we get "n+m+l "n+m m(An Xn-1) X2k + 1 / k - X2k + 1 / k = H x2k + 1 / k 2k+l/k ( 18) In particul a r when n=l, it follows "m+2 ~+1 m "1 AO X2k + 1 / k - X2k + 1 / k = H (X 2k + 1 / k - X2k + 1 / k ) ( 19) Since Hm=O for m~2, then it turns out

"3 A4 A5 X2k + 1 / k = X2k + 1 / k = X2k + 1 / k = .•• and s o for ever. This means that three approximations should be enough to complete the procedure. Now we accept for gra nted the following

A2kP2k/kA2k + B2kMB2k

where

" + / ' P 2k + 1 / k = cov(X 2k + 1- X 2k 1 k "

X2k + 1- X2k + 1 / k

"

"-

P 2k / k = cov(X2k-X2k/k' X2k-X2k/k) M = cov(f - f , f - f ) (which is k k k k obtained from the foregoin g part) It is easy to verify that the solution of Eq. (23) has the form as 1

P

" "'3 X2k+1 / k = X2k + 1 / k

(23)

+ F2kG2F2k

(20)

- >" HiQ(H,)j 2k+l /k-. '-:- 0

(24)

l,J=

where Then the 8q. (14) bec omes after the a pproximation procedure is completed as follows A "

A

A2kX2k / k + B2k U2k + B2kfk

X2k + 1 / k

+ HX" 2k + 1 / k

(21)

It mi g ht be interestin ~ to see that on the right hand side of Eq. (21) there is term HX 2k + 1 / k , but no HX 2k / k as it could be imagined. The reasons for this are following 1. This is a logical consequence of the approximation procedure which is defined above 2. If we replace HX" 2k + 1 / k by that i s to replace b by (~­ 2 in the system of Eqs. (7). If would be the case, the system become unobserva ble one.

"

HX 2k / k '

Y 1 ) w~c

this might

Subtracting Eq. (21) from Eq. (11) a nd taking account of B2kb=B2kfk+HX2k+l' an estimation error equation will be obtained as follows

"

X2k + 1 - X2k + 1 / k A2k (X 2k -

Q = A2kP2k/kA2k+ B2k MB 2k + F2kG2F2k (25)

X2k / k )

+ B2k (f k -

fk)

+ F2kN2k + H(X 2k + 1- '"X2k + 1 / k ) (22) We assume that the quantities in model-l are independent stochastically on that of model-2 as long as we consider these models subjected to various white noises. Then the estimation error covariance matrix can be given as follows (I - H)P2k+l/k(I - H)'

Hitherto, we have obtained the timeupdating equations of the algorithm, that are Eqs. (21) (24). Now, suppose a new measurement has been taken, i.e. Z2k+l' We define an innovation of new measured information as follows (Kailath, 1976)

"

e k + 1 = Z2k+l- Z2k+l/k C2k+l(X2k+l-

X2k + 1 / k )+ Vk + 1

Then the llse of X2k + 1 , given (Z2i) i = 0,1 ,2, ••• ,k+1, can be formulated as follows "

1\

X2k+l/k+l

X2k + 1 / k + P2k+1/kC2k+1

(C2k+1P2k+1/kC2k+1+R2k+l)-1ek+1 (26) The optimality of this estimate can be proved by checking the orthogonality principle, which states that the estimation error must be orthogonal to the innovation process (e ), i=O, i 1,2, ••• Indeed we have

k

E[(X 2k + 1- X2k+1/k+1 )e + 1] = E[(X 2k + 1 - X2k+1/k)ek+J- P 2k + 1 / k C2k + 1 ( C2k+1P2k+1/kC2k+1+ R2k + 1

k

E(e k + 1 e + 1 )

(27)

However, E[(X 2k + 1 - X2k+1/k)ek+1] E{(X 2k + 1-

) -1

X2k + 1 / k ) [(X 2k + 1-

- X2k+1/k)'C2k+1+ V2k + 1]}

Atti t ud e Estimati on o f Earth Observ a ti on SaLe lli tes

(28)

= P2k+1/kC2k+l

k

E(e k + 1 e + 1 ) = E{[C2k+1(X2k+l

- X2K + 1 / k ) -

+ V 2k + 1] [(X 2k + 1

X2k+1/k)'C2k+1

+ V2k + 1 ]}

C2k + 1 P2k+1/k C 2k+1 + R2k + 1 (29) Substituting Eqs. (28) (29) into Eq. (27), we obtain the necessary orthogonality condition fulfilled for the estimate

k J

E [(X 2k + 1 -

X2k+1/k+1 )e + 1

= 0

In summarizing, the recursive estimation algorithm, developed above, will be given as follows A

~ 1 9

tisfactory results of yaw attitude estimation. So, th e conventional kalman filter cannot be applied straightforward. A combined estimation algorithm, proposed in the paper, consists of two parts in according to two variants of system model. The constant disturbance and gyro drift are treated in different ways by these models, and these two models are interconnected through a properly defined approximation procedure in determining the constant disturbance. This approach preserves the observability in both parts, and linear least square estimation algorithm has been given to them. The output of this algorithm is complete estimation of all system variables, including yaw attitude, unknown constant disturbances and gyro drift characteristics. The performance of the algorithm is good enough for many practical needs.

A

A2kX2k/k + B2k U2k + B2kfk

,..

REFERENCES

+ HX 2k + 1 /k i2k+1/k+1 = i2k+1/k + K2k + 1 (Z2k+1 - C2k+1i2k+1/k) P2k+1/kC2k+1(C2k+1P2k+1/kC2k+1

+

R2 Ktl

)-1

1

. . H1Q(H,)J i, j=O

E

P 2k + 1 / k Q

= A2kP2k/kA2k

(30)

+ B2kMB2k + F2kG2F2k

P2k+1/k+1 = P 2k + 1 / k - K2k+1C2k+1P2k+1/k

X20 /_ 1

= X20 :

P 20 /- 1

= P 20

The Eqs. (30) must be complemented by A

Eqs. (13) in order to calculate fk and M. Here the results of a computer simulation of the algorithm are presented in Fig. 1-6. Some typical parameters in the simulation are chosen as follows -5 2 horizon sensor's accuracy:l0 rad sec constant gyro drift:2deg./hr constant disturbance:3 g-cm CONCLUSIONS The yaw attitude estimation of an earth observation satellite is most challenging problem in the case, where earth horizon sensors and rate gyros are being used alone as attitude measuring devices. The constant disrurbance and gyro drift are dominant factors in preventing from achieving sa-

Anderson, B. D. 0., and J. B. Moore (1979). Optimal Filtering. Prentice-Hall, Inc., Englewood Cliffs, N.J. Bowers, J. L., Rodden, J. J., Scott, E. D., and D. B. DeBra (1968). Orbital gyrocomp a ssing heading reference. J. Spacecraft and Rockets,2., 903-910 Bryson, A. E., and W. Kortum (1971). Estimation of the local attitude of orbiting spacecraft. Automatica, 7, no. 2 Kailath, T::- (1 '376) • Lectures on Linear Least-Squares Estimation. Springer Verlag, Wien-New York. Yang, J. C., and C. Q. Suing (1980). Gyrocompassing accuracy improvement by gyro drift estimating in orbit. Paper, presented at The Joint Symposium of Chinese Association of Automation and Chinese Association of Astronautics on Automatic Control in Space, 1980. (in Chinese)

s.

220

H. Yu and Q. P. Chu

2

4

2

o

40

20

Fig. 1.

60

t

(min)

o Fi g . 2.

Roll estimation error covariance in 10- 10 rad 2 vs. time.

40

20

60

t

(min)

Yaw estimation error covariance in 10-8 rad 2 VB. time.

2

o -1

-2

Fig. 3.

Fig.

Estimation error -1 in 10 deg. vs. time in min.

4.

Eetimation error of disturbance b, in 10- 4 g-cm vs. time in min.

4

2

Roll 0

2 t

0 t

-1

-2

-2

-4 Fig . 5.

Estimation error of constant gyro drift in '0-7 rad/sec vs. time.

Fig. 6.

Estimation error of disturbance b 2 in 10- 3 g-cm vs. time in min.