Attitude Determination Using a MEMS-Based Flight Information Measurement Unit

Obtaining precise attitude information is essential for aircraft navigation and control. This paper presents the results of the attitude determination using an in-house designed low-cost MEMS-based flight information measurement unit. This study proposes a quaternion-based extended Kalman filter to integrate the traditional quaternion and gravitational force decomposition methods for attitude determination algorithm. The proposed extended Kalman filter utilizes the evolution of the four elements in the quaternion method for attitude determination as the dynamic model, with the four elements as the states of the filter. The attitude angles obtained from the gravity computations and from the electronic magnetic sensors are regarded as the measurement of the filter. The immeasurable gravity accelerations are deduced from the outputs of the three axes accelerometers, the relative accelerations, and the accelerations due to body rotation. The constraint of the four elements of the quaternion method is treated as a perfect measurement and is integrated into the filter computation. Approximations of the time-varying noise variances of the measured signals are discussed and presented with details through Taylor series expansions. The algorithm is intuitive, easy to implement, and reliable for long-term high dynamic maneuvers. Moreover, a set of flight test data is utilized to demonstrate the success and practicality of the proposed algorithm and the filter design.


Introduction
Obtaining an accurate vehicle attitude is essential for airplane navigation and control applications. The effectiveness of navigation and control is determined by the degree of precision of the navigation and control systems, including inertial measurement units. Traditional units-such as gimbaled gyroscopes, laser gyroscopes, and fiber optic gyroscopes and accelerometers-provide high-precision information for navigational calculations [1]; they are, however, expensive and bulky. With the maturation and advancement of semiconductor manufacturing technology, MEMS sensors are increasingly used in flight attitude calculations [2][3][4][5][6][7][8][9].
The attitude of the aero-vehicle can be determined by integrating the angular rates (pitch, roll, and yaw rates) of the vehicle. Nevertheless, accuracy requirements usually cannot be satisfied by using the inexpensive MEMS sensors. Therefore, some forms of Kalman filtering or complementary fusion algorithms are normally employed to provide more accurate and reliable attitude angles in the MEMS attitude determination systems. The pitch and roll angles can be estimated by using the gravity components from the measurement of the accelerometers [2]. Gravity acceleration cannot be determined accurately from the accelerometer measurements, especially in long-term high dynamic maneuvers. To increase the accuracy and reliability of the attitude determination, a considerable number of advanced fusion technologies using MEMS sensors were proposed in the literature. Using a triad of rate gyros fused with an ultra-short baseline, the differential GPS attitude determination system was presented in [2]. In [4], a gyro-free attitude determination was discussed by using accelerometers, magnetic sensors, and GPS as the primary sensors. A fuzzy logic based closed-loop strapdown attitude system for UAV is presented in [5]. In [7], a constrained Kalman filter was developed to eliminate the effects of kinetic accelerations in the y and z axes in the body frame, while using a differential speed signal to correct the kinetic acceleration of the x axis. A gain-scheduled complementary filter augmented by an acceleration-based switching architecture was proposed in [9] to provide robust performance, even when the vehicle is subject to strong acceleration. A miniature MEMS-based attitude and heading reference system, which makes use of an extended Kalman filter with adaptive gain, and with accelerometers and magnetic sensor as the measurements, was presented in [10]. A precision attitude determination employing a multiple model adaptive estimation scheme was reported in [11]. A Linear Fusion Algorithm for Attitude Determination Using Low Cost MEMS-Based Sensors is considered in [12]. In [13], an integrated MEMS INS/GPS and intelligent artificial neural network with Rauch-Tung-Striebel smoother scheme was proposed for position and orientation determination.
This study considers the utilization of a low-cost MEMS-based flight information measurement unit to determine attitude. The main contribution of this paper is that an attitude estimation algorithm is proposed using MEMs sensors. In the proposed algorithm, an extended Kalman filter is used instead of the classical channel-wise complementary filters for attitude estimation. Specifically, the roll and pitch angles computed from the gravity force decomposition, and the heading angle determined from the electronic compass unit, are selected as the measurement of the Kalman filter. This is not commonly used in airplane attitude determination. The immeasurable gravity accelerations are deduced from the outputs of the three axes accelerometers, the relative accelerations, and the accelerations due to body rotation. The proposed attitude determination algorithm employed the evolution of the four elements of the quaternion method as the dynamical model for the Kalman filter. The constraint of the four elements of the quaternion method is thus treated as a perfect measurement and is integrated into the filter design. By adopting a coupled and time-varying approach, the proposed algorithm performs well under different dynamic maneuvers. The proposed attitude estimation algorithm is intuitive and easy to implement, and due to the vehicle dynamics has been considered in the gravity computation, the algorithm can handle the situation of the long-term high dynamic maneuvers. Detailed derivations establish the approximations of the time varying noise variances of the measured signals.
The proposed algorithm is demonstrated through a set of flight test data collected from the in-house designed MEMS-based attitude determination system. The duration of the flight covers a total of 3,620 s. Details of the evolutions of the attitude angles are examined and explored in the paper. The results show the long-term stability and practicality of the proposed attitude estimation algorithm.

Flight Information Measurement Unit
The determination of attitude considered in this study is based on a set of signals measured from an in house-designed Flight Information Measurement Unit (FIMU) [14] and a GPS receiver. The main function of the FIMU system is to measure the required signals for navigation computation and estimation; this includes the accelerations, angular rates, magnetic fields, differential pressures, static pressure, and temperature in a real-time and continuous manner. Hence, the receiver (an on-board computer) can use this information to calculate the vehicle's attitude, heading, airspeed, altitude, angle of attack, and temperature information. The basic system structure is shown in Figure 1. To avoid or reduce electromagnetic interference during the measurement process (because the magnetic sensor is extremely sensitive to its ambient disturbances), the FIMU is split into two separate units, namely, an inertial and air data measurement unit, and an electronic compass unit. On the one hand, the inertial and air data measurement unit provides the capability for three-axis accelerations, three-axis angular rates, airspeed, altitude, temperature, and angle of attack measurements. The compass system, on the other hand, measures the magnitude of the external magnetic fields. A Microchip PIC micro-controller is used in each measuring unit to perform data acquisition and signal processing, and to transfer the results to the host system for further analysis and processing. An image of the FIMU is shown in Figure 2.

Attitude Determination
Determination of flight attitude involves the computation of aircraft pitch angle, roll angle, and yaw angle. Pitch angle and roll angle can be computed through the measured aircraft accelerations and body rates from the accelerometers and rate gyros. The heading angle is determined by computing the magnetic heading: a magnetic sensor measures the magnetic field, and corrects for the magnetic declination. With the measured three axes acceleration signals and pitch, roll, and yaw rate information -the pitch and roll angles can be determined either by computing the gravitational acceleration components on the body axes, or by using the Euler quaternion method [15]. Computation of the gravitational acceleration components provides long-term accuracy, although it is accompanied by high noise contents. The quaternion method, however, provides low noise contents and fast response to changes in the input signals, but tends to drift with time due to gyro bias errors. The results of pitch and roll computations (derived from gravitational acceleration components and the quaternion method) are mixed by means of complementary filtering, which is performed by implementing sensor data fusion techniques. This sensor data fusion technique can also be incorporated for heading computation by mixing the heading information from both the compass system and from the Euler quaternion method.
For example, assume that at a certain instant, the pitch angle of the aircraft is θ, roll angle is ϕ, and the gravitational force components along the body axes (X-, Y-, and Z-axis) are g x , g y , and g z respectively, as shown in Figure 3 [1]. The relationship between the gravitational acceleration components and the attitude angle are: (1) sin sin cos cos cos Knowing g x , g y , and g z , the roll angle ϕ and pitch angle θ can be computed from: There are alternative ways to obtain ϕ and θ from Equation (1). They are mathematically equivalent. Computation of the gravity vector g g g g is avoided by using Equation (2).
But g x , g y , and g z cannot be measured directly during flight. The accelerations measured from the accelerometers are the total accelerations along the X-, Y-, and Z-axis (a x ,a y ,a z ). The relationship between the measured accelerations (a x ,a y ,a z ) and the gravitational force components, as shown in Figure 4 [1], are: where (U,V,W) are the inertial velocities along the X-, Y-, and Z-axis, whereas (p,q,r) are the pitch rate, roll rate, and yaw rate, respectively. Note that the local level frame is taken as the inertial frame in the study. We may thus derive (g x ,g y ,g z ) by measuring or computing (U,V,W), ( , , ), (p,q,r), and (a x ,a y ,a z ); consequently, these results can be used to determine the aircraft pitch angle θ and roll angle ϕ. The advantage of this method is that no integral operation is involved in the computation. Thus, the measurement error and gyro biases are not accumulated, thus avoiding divergence in attitude computation. It therefore provides the results with long-term accuracy. However, due to the nature of the MEMS accelerometers and gyros and the computation of (U,V,W) and ( , , ), the results usually lead to high noise contents. For magnetic heading computation, assume that the components of the Earth's magnetic field along the X-, Y-, and Z-axis are H X , H Y , and H Z , respectively. Furthermore, the resolved components of H X , H Y , and H Z in the horizontal plane along the heading axis H 1 and at right angles to the heading axis H 2 (as shown in Figure 5) are given by [1]: sin sin cos sin Thus, the magnetic heading of the aircraft ψ M is: If the local magnetic declination is ψ md , the aircraft heading can be determined from ψ = ψ M + ψ md . A common method used to compute the attitude angle is the Euler quaternion method [15], which uses the four symmetrical Euler parameters to define the aircraft attitude. The relationship between the attitude angles and the four parameters are: The four parameters are subjected to the following constraint equation: The dynamics of the four parameters regarding the aircraft body rate (p,q,r) are characterized in the following form: (9) Therefore, the attitude angle in Equation (7) may be calculated by solving Equation (9). Only body information-(p,q,r) directly measured from rate gyros-is required to solve Equation (9). However, when MEMS gyros are used for body rate measurement, long-term drift is usually encountered due to gyro bias errors and integral operation for solving Equation (9).

Kalman Filter Design
As discussed in the previous section, the attitude computed from the gravitational force components provides long-term accuracy with high noise contents. The Euler quaternion method, however, provides less noisy results but suffers from a long-term drift problem. Either method alone may be inadequate for attitude computation. Hence, in this study, the extended Kalman filter [16] is implemented to integrate the attitude computation from the gravitational force components and from the Euler quaternion method.
The Kalman filter is a model-based estimation technique. The dynamics that characterize the relationship between the aircraft body rate (pitch, roll, and yaw rates) and the four parameters of the quaternion method described in Equation (9) comprises the dynamic model of the Kalman filter. The relationship between the attitude angle and the four parameters in Equation (7) is chosen as the output equation of the filter. Thus, the dynamic system can be expressed as: T is the state of the system, y = [ϕ θ ψ 1] T is the output, and w and v are the precess and measurement noise, respectively. The system matrix A and output function h(x) are defined as: The measurements of pitch angle θ and roll angle ϕ are determined from the gravitation force components Equation (1), while the heading angle ψ is determined from the measured magnetic heading Equation (6), and is corrected with local magnetic declination. The gravitation force components (g x ,g y ,g z ) are determined from Equation (3), that is: Except the measurement of accelerations (a x ,a y ,a z ) and body rates (p,q,r), the velocities (U,V,W) and velocities' rates ( , , ) along the body axes are also required. This information can be determined from: (13) where V T is the airspeed, α is the angle of attack, and α is the side slip angle. This approach requires an angle of attack sensor and a sideslip sensor to acquire the angle of attack and angle of sideslip information. Equation (13) is valid only for no wind is present. Therefore, in this study, the inertial velocity components (U,V,W) are determined from the ground speed information (V N , the north velocity, V E , the east velocity, and V D , the down velocity of the aircraft) from the GPS receiver through yaw, pitch, and roll rotation sequence. Equation (14) expresses the coordinate transformation for (V N , V E , V D ) and (U,V,W): (1 tan tan ) (1 tan tan ) tan sin sin cos sin cos cos sin sin sin cos cos cos The velocities' rates ( , , ) are computed from ( , , ) through the coordinate transformation as in Equation (14) for (U,V,W). The velocities' rates ( , , ) are computed from ( , , ) by the linear approximation 1 ∆ , 1 ∆ , and 1 ∆ . The angle of attack and sideslip sensors are not installed on the airplane when conducting the flight test for this research. Therefore, the results of using Equation (13) to compute the inertial velocity components are not investigated in this paper.
For real time computation, the dynamical system in Equation (10) is expressed in discrete time representation as: (15) where ∆ , which can be represented by the following form [17]: with ΔP(n) = p(n)Δt, ΔQ(n) = q(n)Δt, ΔR(n) = r(n)Δt, Δt is the sampling time (0.05 s in this study), p(n), q(n), r(n) are the pitch rate, roll rate, and yaw rate at time index n. The output functionh[x(n)]-is nonlinear. Therefore, the extended Kalman filter is chosen to perform the filtering for the attitude determination in this study. The computation process is depicted in Figure 6, and the Kalman filtering details are shown in Figure 7. In the estimation process, after initialization, the iteration loop starts with computing the process matrix F(n−1) and the process noise variance, after which the time updates of the state estimate and estimation-error covariance are performed. Subsequently, the linearized measurement function H n and the variances of the measurement noises are determined. The function H n is obtained by taking the derivatives of the observation function h(x) with respect to the state of the filter. The state of the filter is the quaternion components. Therefore, the function H n is given in the following form: (17) [ ] With all the necessary information available, the last step of the iteration loop is to perform the measurement update of the state estimate and estimation-error covariance . Details of the process and measurement noises are presented in the following section.

Noise Characteristics
In the Kalman filter design, we use the covariances of the process noise and measurement noise as design parameters to achieve minimum variance of the estimation error. Hence, we must determine the variances of the process and measurement noises as precisely as possible. For the attitude determination considered in this study, the process noises are mainly derived from the outputs of the gyros (p,q,r). Assuming that that , , and ̃, with , , and the mean of p, q, and r respectively, and , , and ̃ are the deviations of p, q, and r, Equation (9) can be written as: The second term on the right hand side is considered as the process noise w(t) in this study. In the sequel, represents the expectation or mean of the signal , and indicates the deviation from the mean value.
For a discrete time system, the process noise w(n) is: It is simple to show that the variance of the process noise is with: where , , and are the variances of p, q, and r, respectively. The characteristics of the sources of the measurement determine the variances of the measurement noises. In this study, the pitch angle θ and roll angle ϕ are computed from the measured gravity components in Equations (1) and (12), while the heading angle is determined from the electronic compass system in Equation (6). We can therefore represent the variance of the measurement noises in the following form: represents the variance of pitch and roll angles (θ,ϕ); is the variance of heading angle; is the variance of the constraint in Equation (8), and is equal to zero. To avoid numerical problem, the variance can be set to a small constant instead of zero. Since the pitch and roll angles are computed from Equation (12), the deviations of g x , g y , and g z , denoted as g , g , and g , can be decided by expanding Equation (12) in a Taylor series about the mean of the measured variables, (a x ,a y ,a z ), (U,V,W), ( , , ), (p,q,r) and neglecting the second order terms to obtain: where: The roll angle ϕ and pitch angle θ are computed from: Expanding the function in the Taylor series about the mean values and neglecting the second and higher order terms, the deviation of can be approximated by: To determine , the function tan is also expanded in a Taylor series, and the deviation is approximated as: Similarly, the deviation of can be approximated by: g cos g cos g g g cos g g g sin g cos g g g cos g g g sin g 0 g g g g g g g g g cos g g sin g g g cos g g g sin g g g g g g (27) x y x y z We also assume that , , and are all uncorrelated. The variance of the heading angle is equal to since the only difference between the heading angle ψ and the magnetic heading ψ M is the magnetic declination angle.

Design Evaluation
To evaluate the design of the flight information measurement unit and the extended Kalman filtering for attitude determination, the system is tested on an ultralight aircraft. The test dataincluding the three axes body rates (p,q,r), accelerations (a x ,a y ,a z ), magnetic field (H X ,H Y ,H Z ), and information from GPS receiver-are recorded during flight test. The flight was performed in Central Taiwan at approximately 120.44 degrees East Longitude and 23.77 degrees North Latitude area. The flight trajectory is shown in Figure 8. The evolutions of the attitude, including pitch, roll, and heading angles are shown in Figure 9. The attitude angles (directly generated from the measured data and computed after the Kalman filtering) are both shown in Figure 9. The duration of the flight covers a total of 3,620 s. The results are further examined and discussed below. In particular, during certain time intervals, the aircraft performed specific flight trajectories: during the time interval 1,106-1,130 s, the aircraft was making a left turn; from 1,420-1,430 s, the aircraft ascended; the aircraft descended during the time interval 2,340-2,380 s; during the time interval 2,600-2,675 s, the aircraft was making a right turn; from 2,745-2,755 s, the aircraft was in level flight with wing wiggling; and from 2,810-2,900 s, the aircraft was in S-turn maneuver. Much of the noisy signals of the attitude angles in the followings figures are the results deduced from the measurement, whereas the less noisy signals are the results of the filtered signals.  The results for nearly level flight with wing wiggling (intentionally performed) are shown in Figure 10. In this phase, the altitude was sustained at approximately 392 to 393 meters height, as shown in Figure 10(a). The roll angle, as shown in Figure 10(c), responded as commanded. Figure 10(d) shows the tendency of turning to the left. The average of the roll angles in Figure 10(c) is negative, which confirms that the airplane is turning left slowly. The results also show that the time delay of the attitude estimation is approximately 0.5 s due to filter computation. Figure 11 indicates the results for the climbing up condition. During the time period of 1,420 to 1,430 s, the airplane climbed from 1,009 m to 1,019 m height. The results during the descent flight (time duration 2,340 to 2,380 s) are shown in Figure 12. In this descent flight phase, the altitude decreases from 290 to 190 m height. Evolutions of the pitch angle in Figure 12(b) are all positive, with an average of roughly 10 degrees (which is normal for aircraft descent flight). Figure 13 shows the results when the airplane made a left turn. The data from the GPS receiver, as shown in Figure 13(a), indicate that the airplane is turning left. When making a left turn, the airplane banks to the left before starting to turn. The results in Figure 13(c), negative roll angles, show that the airplane rolls to the left, which confirms that the airplane is turning to the left. The heading angle changes from approximately 200 degrees to 60 degrees in Figure 13(d), which verifies that the airplane is making a left turn. Figure 14 presents the results when the airplane made a right turn. In the right turning phase, the airplane has a positive roll angle, which is evident in Figure 14(c). During the time period from 2,810 to 2,900 s, the airplane made an S-turn. Evolutions of the attitude angles during this maneuver are shown in Figure 15. Figure 15(a) shows the trajectory of the S-turn maneuver from the GPS data. During this particular maneuver, the airplane makes a left turn before following with a right turn. Figure 15(c) shows that the roll angle turns to a negative value first, then returns to zero degrees for the left maneuver. Thereafter, it turns positive before returning to zero degrees for a right turning operation. The heading angle shown in Figure 15(d) starts from approximately 340 degrees and gradually decreases to roughly 25 degrees, then increases back to 335 degrees, which correctly corresponds to the S-turn operation.   From the above examinations, we demonstrated that the proposed attitude determination algorithm estimates the attitude angles successfully. The algorithm is also reliable for long-term and high-dynamic maneuvers.
To further investigate the performance of the proposed attitude estimation method, the innovation around ±3 standard deviations (extracted from the innovation covariance) and the quaternion norm constraint in Equation (8) are examined. Figure 16(a-c) represent the results of the innovation when the airplane made a left turn. Figure 16(d-f) are the results for nearly level flight with wing wiggling. At some points, the innovation exceeds the boundaries of the ±3 standard deviations. This is due to the fact that the extended Kalman filter is an approximate solution. It relies on the linearization to propagate the mean and covariance of the state. The norms of the quaternion computed from the extended Kalman filter for all over the duration of flight are displayed in Figure 17. The results show that the constraint optimization is well performed.
Because there is no primary standard, the proposed attitude estimation algorithm is also tested using a motion platform under a controlled environment. An image of the test arrangement is shown in Figure 18. A cantilever beam is mounted atop the motion platform. A specially designed movable mounting plate is constructed at the end of the beam. The IMU is placed on plate. The mounting plate is arranged in such a way that the IMU is mounted with 0 degree pitch angle and −20 degrees roll angles. An on-board computer is mounted on top of the center of the motion platform to record the outputs from the IMU during the test.   accelerometers and the body rates from the gyro outputs, the pitch angle and roll angle are computed by using each of the gravity force components, the Euler quaternion method (which can be considered as a direct 6-DOF mechanization for attitude determination [18]), and the proposed extended Kalman filtering. The results are shown in Figure 19. The results in Figure 19 indicate that the gravity force decomposition will provide accurate results of 0 degree pitch angle and −20 degrees roll angle but suffer from noisy contents. The roll angle is correctly deduced from the quaternion computation with less noisy contents. However, the pitch angle tends to diverge. The results using the proposed extended Kalman filtering provide both accurate and less noisy results for the pitch and roll angles.

Conclusions
This study presented attitude determination using a low cost MEMS-based flight information measurement unit. The proposed attitude estimation algorithm is intuitive, easy to implement, and reliable for long-term high dynamic maneuvers. The algorithm utilized the evolution of the four elements of the quaternion method as the dynamical model for the Kalman filter. The roll and pitch angles (computed from the gravity force decomposition) and the heading angle (determined from the electronic compass unit) are selected as the measurements of the Kalman filter. The immeasurable gravity accelerations are deduced from the outputs of the three axes accelerometers, the relative accelerations, and the accelerations due to body rotation. The constraint of the four elements of the quaternion method is treated as perfect measurement and is integrated into the filter design. Approximations of the time varying noise variances of the measured signals are established with detailed derivations. The proposed algorithm is successfully demonstrated through a set of flight test data collected from the in house-designed MEMS based attitude determination system. To improve the performance of the attitude estimation, estimation of the gyro drift and velocity derivatives can be included into the Kalman filter formulation with the expense of computation demands.