# Unscented Kalman Filter Applied to the Spacecraft Attitude Estimation with Euler Angles

Authors: Roberta Veloso Garcia; Helio Koiti Kuga; Maria Cecilia F. P. S. Zanardi
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The aim of this work is to test an algorithm to estimate, in real time, the attitude of an artificial satellite using real data supplied by attitude sensors that are on board of the CBERS-2 satellite (China Brazil Earth Resources Satellite). The real-time estimator used in this work for attitude determination is the Unscented Kalman Filter. This filter is a new alternative to the extended Kalman filter usually applied to the estimation and control problems of attitude and orbit. This algorithm is capable of carrying out estimation of the states of nonlinear systems, without the necessity of linearization of the nonlinear functions present in the model. This estimation is possible due to a transformation that generates a set of vectors that, suffering a nonlinear transformation, preserves the same mean and covariance of the random variables before the transformation. The performance will be evaluated and analyzed through the comparison between the Unscented Kalman filter and the extended Kalman filter results, by using real onboard data.


---

## Body

## 1. Introduction

The attitude of a spacecraft is defined by its orientation in space related to some reference system. The importance of determining the attitude is related not only to the performance of attitude control but also to the precise usage of information obtained by payload experiments performed by the satellite. The attitude estimation is the process of calculating the orientation of the spacecraft in relation to a reference system from data supplied by attitude sensors. Chosen the vector of reference, an attitude sensor measures the orientation of these vectors with respect to the satellite reference system [1]. Once these one or more vectors measurements are known, it is possible to compute the orientation of the satellite processing these vectors, using methods of attitude estimation. There are several methods for determining the attitude of a satellite. Each method is appropriate to a particular type of application and meets the needs such as available time for processing and accuracy to be attained. However, all methods need observations that are obtained by means of sensors installed on the satellite. The sensors are essential for attitude estimation, because they measure its orientation relative to some referential, for example, the Earth, the sun, or a star. In this work, the satellite attitude is described by Euler angles, due to its easy geometric interpretation, and the method to estimate the attitude used is the Unscented Kalman Filter. This method is capable of performing state estimation in nonlinear  systems, besides taking into account measurements provided by different attitude sensors. In this work there were considered real data supplied by gyroscopes, infrared Earth sensors, and digital sun sensors. These sensors are on board of the CBERS-2 satellite (China-Brazil Earth Resources Satellite), and the measurements were collected by the Satellite Control Centre of INPE (Brazilian Institute for Space Research).


## 2. Representation of Attitude by Euler Angles

The attitude of an artificial satellite is directly related to its orientation in space. Through the attitude one can know the spatial orientation of the satellite, since in most cases it can be considered as a rigid body, where the attitude is expressed by the relationship between two coordinate systems, one fixed on the satellite and another associated with a reference system, for example, inertial system [2]. For a good performance of a mission it is essential that the satellite be stabilized in relation to a specified attitude. The attitude stabilization is achieved by the on-board attitude control, which is designed to acquire and maintain the satellite in a predefined attitude. The CBERS-2 attitude is stabilized in three axes nominally geo-pointed and can be described with respect to the orbital system. In this reference system, the movement around the direction of the orbital velocity is called roll. The movement around the direction normal to the orbit is called pitch, and finally the movement around the direction Nadir/Zenith is called yaw. To transform a vector represented in a given reference to another it is necessary to define a matrix of direction cosines (R), where its elements are written in terms of Euler angles (ϕ, θ, ψ) [3]. The rotation sequence used in this work for the Euler angles was the 3-2-1, where the coordinate system fixed in the body of the satellite (x, y, z) is related to the orbital coordinate system (xo, yo, zo) through the following sequence of rotations: (i)
1st rotation of an angleψ (yaw angle) around the zo axis,(ii)
2nd rotation of an angleθ (roll angle) around an intermediate axis y′,(iii)
3rd rotation of an angleϕ (pitch angle) around the x-axis.The matrix obtained through the 3-2-1 rotation sequence is given by(2.1)R=[C(θ)C(ψ)C(θ)S(ψ)-S(θ)S(ϕ)S(θ)C(ψ)-S(ψ)C(ϕ)S(ϕ)S(θ)S(ψ)+C(ϕ)C(ψ)S(ϕ)C(θ)S(θ)C(ψ)C(ϕ)+S(ϕ)S(ψ)C(ϕ)S(θ)S(ψ)-S(ϕ)C(ψ)C(θ)C(ϕ)],
where R is the matrix of direction cosines with S=sin, C=cos, and T=tan.By representing the attitude of a satellite with Euler angles, the set of kinematic equations are given by [4](2.2)[ϕ̇θ̇ψ̇]=[1S(ϕ)T(θ)C(ϕ)T(θ)0C(ϕ)-S(ϕ)0S(ϕ)C(θ)C(ϕ)C(θ)]{[ŵxŵyŵz]-R[0-ω00]},
where ω0 is the orbital angular velocity and ŵx, ŵy, and ŵz are the components of the angular velocity on the satellite system.


## 3. The Measurements System of Satellite

In order to estimate the satellite attitude accurately, several types of sensors, including gyros, earth sensors, and solar sensors, are used in the attitude determination system. The equations of these sensors are introduced here.


### 3.1. The Model for Gyros

The advantage of a gyro is that it can provide the angular displacement and/or angular velocity of the satellite directly. However, gyros have an error due to drifting, meaning that their measurement error increases with time. In this work, the rate-integration gyros (RIGs) are used to measure the angular velocities of the roll, pitch, and yaw of the satellite. The mathematical model of the RIGs is [4](3.1)ΔΘi=∫0Δt(ωi+εi)dt(i=x,y,z),
where ΔΘ are the angular displacement of the satellite in a time interval Δt, and εi are components of bias of the gyroscope.Thus, the measured components of the angular velocity of the satellite are given by(3.2)ω⃗̂=(ΔΘ⃗Δt)-ε⃗̂-η⃗1=g⃗-ε⃗-η⃗att,
where g⃗(t) is the output vector of the gyroscope, and η⃗1(t) represents a Gaussian white noise process covering all the remaining unmodelled effects:(3.3)E[ηattηattT]=σatt2.


### 3.2. The Measurement Model for Infrared Earth Sensors (IRESs)

One way to compensate for the drifting errors present in gyros is to use the earth sensors. These sensors are located on the satellite and aligned with their axes ofroll and pitch. In the work, two earth sensors are used, with one measuring the roll angle and the other measuring the pitch angle. In principle, an earth sensor cannot measure the yaw angle. The measurement equations for the earth sensors are given as [5](3.4)ϕH=ϕ+ηϕH,θH=θ+ηθH,
where ηϕH and ηθH are the white noise representing the small remaining misalignment, installation, and/or assembly errors assumed to be gaussian:(3.5)E[ηϕHηϕHT]=E[ηθHηθHT]=σIRES2.


### 3.3. The Measurement Model for Digital Solar Sensors (DSSs)

Since an earth sensor is not able to measure theyaw angle, the solar sensors are used by the Attitude Control System in order to overcome this problem. However, these sensors do not provide direct measurements but coupled angle of pitch (αθ) and yaw (αψ). The measurement equations for the solar sensor are obtained as follows [5]:(3.6)αψ=tan-1(-SySxcos⁡(60°)+Szcos⁡(150°))+ηαψ
when |Sxcos(60°)+Szcos(150°)|≥cos(60°), and(3.7)αθ=24°-tan-1(SxSz)+ηαθ
when |24°-tan-1(Sx/Sz)|<60°, where ηαψ and ηαθ are the white noise representing the small remaining misalignment, installation, and/or assembly errors assumed Gaussian:(3.8)E[ηαψηαψT]=E[ηαθηαθT]=σDSS2.
The conditions are such that the solar vector is in the field of view of the sensor, and Sx, Sy, and Sz are the components of the unit vector associated to the sun vector in the satellite system and given by(3.9)Sx=S0x+ψ̂S0y-θ̂S0z,Sy=S0y-ψ̂S0x+ϕ̂S0z,Sz=S0z-ϕ̂S0y+θ̂S0z,
where S0x, S0y, and S0z are the components of the sun vector in the orbital coordinate system and ϕ̂, θ̂, and ψ̂ are the Euler angles-estimated attitude.


## 3.1. The Model for Gyros

The advantage of a gyro is that it can provide the angular displacement and/or angular velocity of the satellite directly. However, gyros have an error due to drifting, meaning that their measurement error increases with time. In this work, the rate-integration gyros (RIGs) are used to measure the angular velocities of the roll, pitch, and yaw of the satellite. The mathematical model of the RIGs is [4](3.1)ΔΘi=∫0Δt(ωi+εi)dt(i=x,y,z),
where ΔΘ are the angular displacement of the satellite in a time interval Δt, and εi are components of bias of the gyroscope.Thus, the measured components of the angular velocity of the satellite are given by(3.2)ω⃗̂=(ΔΘ⃗Δt)-ε⃗̂-η⃗1=g⃗-ε⃗-η⃗att,
where g⃗(t) is the output vector of the gyroscope, and η⃗1(t) represents a Gaussian white noise process covering all the remaining unmodelled effects:(3.3)E[ηattηattT]=σatt2.


## 3.2. The Measurement Model for Infrared Earth Sensors (IRESs)

One way to compensate for the drifting errors present in gyros is to use the earth sensors. These sensors are located on the satellite and aligned with their axes ofroll and pitch. In the work, two earth sensors are used, with one measuring the roll angle and the other measuring the pitch angle. In principle, an earth sensor cannot measure the yaw angle. The measurement equations for the earth sensors are given as [5](3.4)ϕH=ϕ+ηϕH,θH=θ+ηθH,
where ηϕH and ηθH are the white noise representing the small remaining misalignment, installation, and/or assembly errors assumed to be gaussian:(3.5)E[ηϕHηϕHT]=E[ηθHηθHT]=σIRES2.


## 3.3. The Measurement Model for Digital Solar Sensors (DSSs)

Since an earth sensor is not able to measure theyaw angle, the solar sensors are used by the Attitude Control System in order to overcome this problem. However, these sensors do not provide direct measurements but coupled angle of pitch (αθ) and yaw (αψ). The measurement equations for the solar sensor are obtained as follows [5]:(3.6)αψ=tan-1(-SySxcos⁡(60°)+Szcos⁡(150°))+ηαψ
when |Sxcos(60°)+Szcos(150°)|≥cos(60°), and(3.7)αθ=24°-tan-1(SxSz)+ηαθ
when |24°-tan-1(Sx/Sz)|<60°, where ηαψ and ηαθ are the white noise representing the small remaining misalignment, installation, and/or assembly errors assumed Gaussian:(3.8)E[ηαψηαψT]=E[ηαθηαθT]=σDSS2.
The conditions are such that the solar vector is in the field of view of the sensor, and Sx, Sy, and Sz are the components of the unit vector associated to the sun vector in the satellite system and given by(3.9)Sx=S0x+ψ̂S0y-θ̂S0z,Sy=S0y-ψ̂S0x+ϕ̂S0z,Sz=S0z-ϕ̂S0y+θ̂S0z,
where S0x, S0y, and S0z are the components of the sun vector in the orbital coordinate system and ϕ̂, θ̂, and ψ̂ are the Euler angles-estimated attitude.


## 4. Attitude Estimation Methods

The goal of an estimator is to calculate the state vector (attitude) based on a set of observations (sensors) [6]. In other words, it is an algorithm capable of processing measurements to produce, according with a given criterium, a minimum error estimate of the state of a system. In this paper, the real-time estimator used to estimate the satellite attitude is a variant of the Kalman filter, applied to problems that present some nonlinearity. This estimator is described as follows.


### 4.1. Unscented Kalman Filter

The basic premise behind the Unscented Kalman Filter (UKF) is that it is easier to approximate a Gaussian distribution than it is to approximate an arbitrary nonlinear function. Instead of linearizing to first order using Jacobian matrices, the UKF uses a rational deterministic sampling approach to capture the mean and covariance estimates with a minimal set of sample points. The nonlinear function is applied to each point, in turn, to yield a cloud of transformed points. The statistics of the transformed points can then be calculated to form an estimate of the nonlinearly transformed mean and covariance. We present an algorithmic description of the UKF omitting some theoretical considerations, left to [7, 8].Consider the system model given by(4.1)ẋk+1=f(xk,k)+Gkwk,yk=h(xk,k)+νk,
where xk is the n×1 state vector and yk is the m×1 measurement vector. We assume that the process noise wk and measurement-error noise νk are zero-mean Gaussian noise processes with covariances given by Qk and Rk, respectively. In this work the state vector at time k is defined by the Euler angles and gyro biases:(4.2)x̂k=[ϕ,θ,ψ,εx,εy,εz]T.Performing the necessary simplifications (small Euler angles) in the set of (2.2), the attitude angles and gyro angular velocity biases are modelled as follows:(4.3)ϕ̇(t)=ω0sin⁡ψ̂+ω̂x+θ̂ω̂z,θ̇(t)=ω0cos⁡ψ̂+ω̂y+ϕ̂ω̂z,ψ̇(t)=ω0(θ̂sin⁡ψ̂-ϕ̂cos⁡ψ̂)+ω̂z+ϕ̂ω̂y,ε⃗̇(t)=0.Given the state vector at stepk-1, we compute a collection of sigma-points, stored in the columns of the n×(2n+1) sigma point matrix χk-1, where n is the dimension of the state vector. In our case, n=6, so χk-1 is a 6×13 matrix. The columns of χk-1 are computed by(4.4)(χk-1)0=x̂k-1,(χk-1)i=x̂k-1+((n+λ)Pk-1)i,i=1,…,n,(χk-1)i=x̂k-1-((n+λ)Pk-1)i-n,i=n+1,…,2n,
in which λϵℜ, ((n+λ)Pk-1)i is the ith column of the matrix square root of (n+λ)Pk-1.Onceχk-1 computed, we perform the prediction step by first predicting each column of χk-1 through time by Δt using(4.5)(χk)i=f((χk-1)i),i=0,…,2n,
where f is differential equation defined in (2.2) or (4.3). In our formulation, we perform a numerical Runge-Kutta integration.With(χk)i being calculated, the a priori state estimate is(4.6)x̂k-=∑i=02nWim(χk)i,
where Wim are weights defined by(4.7)W0m=λ(n+λ),Wim=12(n+λ)i=1,…,2n.As the last part of the prediction step, we calculate thea priori error covariance with(4.8)Pk-=∑i=02nWic[(χk)i-x̂k-][(χk)i-x̂k-]T+Qk,
where Qk is the process error covariance matrix, and the weights are defined by(4.9)W0c=λ(n+λ)+(1-α2+β2),Wic=12(n+λ)i=1,…,2n,
where α is a scaling parameter which determines the spread of the sigma points and β is a parameter used to incorporate any prior knowledge about the distribution of x [9].To compute the correction step, we first must transform the columns ofχk through the measurement function to Yk. In this way(4.10)(Yk)i=h((χk)i),i=0,…,2n,ŷk-=∑i=02nWim(Yk)i.With the mean measurement vectorŷk-, we compute the a posteriori state estimate using(4.11)x̂k=x̂k-+Kk(yk-ŷk-),
where Kk is the Kalman gain. In the UKF formulation, Kk is defined by(4.12)Kk=Px̂k,ŷkPŷk,ŷk-1,
where(4.13)Pŷk,ŷk=∑i=02nWi(c)[(Yk)i-ŷk-][(Yk)i-ŷk-]T+Rk,Px̂k,ŷk=∑i=02nWi(c)[(χk)i-x̂k-][(Yk)i-ŷk-]T,
where Rk represents the measurement error covariance matrix.Finally, the last calculation in the correction step is to compute thea posteriori estimate of the error covariance given by(4.14)Pk=Pk--KkPŷk,ŷkKkT.


## 4.1. Unscented Kalman Filter

The basic premise behind the Unscented Kalman Filter (UKF) is that it is easier to approximate a Gaussian distribution than it is to approximate an arbitrary nonlinear function. Instead of linearizing to first order using Jacobian matrices, the UKF uses a rational deterministic sampling approach to capture the mean and covariance estimates with a minimal set of sample points. The nonlinear function is applied to each point, in turn, to yield a cloud of transformed points. The statistics of the transformed points can then be calculated to form an estimate of the nonlinearly transformed mean and covariance. We present an algorithmic description of the UKF omitting some theoretical considerations, left to [7, 8].Consider the system model given by(4.1)ẋk+1=f(xk,k)+Gkwk,yk=h(xk,k)+νk,
where xk is the n×1 state vector and yk is the m×1 measurement vector. We assume that the process noise wk and measurement-error noise νk are zero-mean Gaussian noise processes with covariances given by Qk and Rk, respectively. In this work the state vector at time k is defined by the Euler angles and gyro biases:(4.2)x̂k=[ϕ,θ,ψ,εx,εy,εz]T.Performing the necessary simplifications (small Euler angles) in the set of (2.2), the attitude angles and gyro angular velocity biases are modelled as follows:(4.3)ϕ̇(t)=ω0sin⁡ψ̂+ω̂x+θ̂ω̂z,θ̇(t)=ω0cos⁡ψ̂+ω̂y+ϕ̂ω̂z,ψ̇(t)=ω0(θ̂sin⁡ψ̂-ϕ̂cos⁡ψ̂)+ω̂z+ϕ̂ω̂y,ε⃗̇(t)=0.Given the state vector at stepk-1, we compute a collection of sigma-points, stored in the columns of the n×(2n+1) sigma point matrix χk-1, where n is the dimension of the state vector. In our case, n=6, so χk-1 is a 6×13 matrix. The columns of χk-1 are computed by(4.4)(χk-1)0=x̂k-1,(χk-1)i=x̂k-1+((n+λ)Pk-1)i,i=1,…,n,(χk-1)i=x̂k-1-((n+λ)Pk-1)i-n,i=n+1,…,2n,
in which λϵℜ, ((n+λ)Pk-1)i is the ith column of the matrix square root of (n+λ)Pk-1.Onceχk-1 computed, we perform the prediction step by first predicting each column of χk-1 through time by Δt using(4.5)(χk)i=f((χk-1)i),i=0,…,2n,
where f is differential equation defined in (2.2) or (4.3). In our formulation, we perform a numerical Runge-Kutta integration.With(χk)i being calculated, the a priori state estimate is(4.6)x̂k-=∑i=02nWim(χk)i,
where Wim are weights defined by(4.7)W0m=λ(n+λ),Wim=12(n+λ)i=1,…,2n.As the last part of the prediction step, we calculate thea priori error covariance with(4.8)Pk-=∑i=02nWic[(χk)i-x̂k-][(χk)i-x̂k-]T+Qk,
where Qk is the process error covariance matrix, and the weights are defined by(4.9)W0c=λ(n+λ)+(1-α2+β2),Wic=12(n+λ)i=1,…,2n,
where α is a scaling parameter which determines the spread of the sigma points and β is a parameter used to incorporate any prior knowledge about the distribution of x [9].To compute the correction step, we first must transform the columns ofχk through the measurement function to Yk. In this way(4.10)(Yk)i=h((χk)i),i=0,…,2n,ŷk-=∑i=02nWim(Yk)i.With the mean measurement vectorŷk-, we compute the a posteriori state estimate using(4.11)x̂k=x̂k-+Kk(yk-ŷk-),
where Kk is the Kalman gain. In the UKF formulation, Kk is defined by(4.12)Kk=Px̂k,ŷkPŷk,ŷk-1,
where(4.13)Pŷk,ŷk=∑i=02nWi(c)[(Yk)i-ŷk-][(Yk)i-ŷk-]T+Rk,Px̂k,ŷk=∑i=02nWi(c)[(χk)i-x̂k-][(Yk)i-ŷk-]T,
where Rk represents the measurement error covariance matrix.Finally, the last calculation in the correction step is to compute thea posteriori estimate of the error covariance given by(4.14)Pk=Pk--KkPŷk,ŷkKkT.


## 5. Results

Here, the results and the analysis from the algorithm developed to estimate the attitude are presented. To validate and to analyze the performance of the estimators, real sensors data from the CBERS-2 satellite were used (see [10, 11]). The CBERS-2 satellite was launched on October 21, 2003. The measurements are for the month of April 2006, available to the ground system at a sampling rate of about 8.56 sec. The algorithm was implemented through MATLAB software. To check the performance the UKF, their results were compared with the estimated attitude by the more conventional EKF (Extended Kalman Filter), considering the following set of initial conditions:(i)
initial attitude:ϕ=θ=ψ=0 deg;(ii)
initialbias of gyro: εx=5.56 deg/hour, and εy=0.87 deg/hour, εz=6.12 deg/hour;(iii)
initial covariance (P): σϕ,θ,ψ2 = (0.5 deg)2 (error related to the attitude), and σbg2 = (1 deg/hour)2 (error related to the drift of gyro);(iv)
observation noise covariance (R): σDSS2 = (0.3 deg)2 (sun sensor), are σIRES2 = (0.03 deg)2 (earth sensor);(v)
dynamic noise covariance (Q): σatt2 = (0.1 deg)2 (noise related to the attitude), σDbgx2=σDbgy2 = (0.01 deg/hour)2, and σDbgz2 = (0.005 deg/hour)2 (noise related to the drifting of gyro).The real measurements obtained by the attitude sensors (digital sun sensors, infrared Earth sensors, and gyroscopes) are shown in Figure1.Real measurements supplied by attitude sensors.
(a)(b)(c)In Figures2 and 3 are observed the behavior of attitude and the biases of gyro during the period analyzed. The average estimated values for the axes of roll and pitch, considering the Unscented Kalman Filter, are in the order of −0.47 deg and −0.45 deg, respectively, and their standard deviations are about 0.02 deg. For the yaw axis the estimate seems not to behave randomly and its average estimated value is about −1.47 deg with standard deviation 0.3 deg. The attitude estimated by the extended Kalman filter had their values for the axes roll and pitch in the order of about −0.49 deg and −0.43 deg with standard deviation about 0.05 deg. For the axis pitch, its average value is −1.42 deg and standard deviation is 0.36 deg.Attitude estimated by the Unscented Kalman filter end Extended Kalman filter.
(a)(b)(c)Bias estimated by the Unscented Kalman filter end Extended Kalman filter.
(a)(b)(c)Figures4 and 5 present the standard deviations for both estimators for the attitude and the bias of the gyro. It is observed that the attitude standard deviations and the standard deviations of the gyro bias decrease with a tendency to stabilize around a value for both filters. However, the graphs show the superiority of UKF, because in most cases it works within a range of protection better than the EKF.Attitude errors (Standard Deviation) estimated.
(a)(b)(c)Bias errors (Standard Deviation) estimated.
(a)(b)(c)In Figures6 and 7, we can see that the residues of sun sensors and Earth sensors have the same behavior for both estimators. For the Earth sensors, the residuals obtained by the both estimators are smaller and show a tendency to zero mean. However, the residues of UKF are still lower, being at about −0.009 deg for IRES1 and 0.004 deg for IRES2. Already the residues of the EKF are approximately 0.01 deg and −0.027 deg for IRES1 and IRES2, respectively.Residuals related to DSS attitude sensors.
(a)(b)Residuals related to IRES attitude sensors.
(a)(b)These results seem consistent because in this case it is not possible to compare the estimated values with true values, since these values are not known.


## 6. Final Comments

In this paper, the Unscented Kalman Filter estimator applied in nonlinear systems was presented for use in real-time attitude estimation. The main objective was to estimate the attitude of a CBERS-2 like satellite, using real data provided by sensors that are on board of the satellite. To verify the consistency of the estimator, the attitude was estimated by two different methods. The usage of real data from on-board attitude sensors poses difficulties like mismodelling, mismatch of sizes, misalignments, unforeseen systematic errors, and postlaunch calibration errors. However, it is observed that, although the EKF and UKF have roughly the same accuracy, the UKF leads to a convergence of the state vector faster than the EKF. This fact was expected, since the UKF prevents the linearizations needed for EKF, when the system has some nonlinearity in their equations.



---

*Source: 985429-2011-12-14.xml*