Effect of Strapdown Integration Order and Sampling Rate on IMU-Based Attitude Estimation Accuracy

This paper deals with the strapdown integration of attitude estimation Kalman filter (KF) based on inertial measurement unit (IMU) signals. In many low-cost wearable IMU applications, a first-order is selected for strapdown integration, which may degrade attitude estimation performance in high-speed angular motions. The purpose of this research is to provide insights into the effect of the strapdown integration order and sampling rate on the attitude estimation accuracy for low-cost IMU applications. Experimental results showed that the effect of integration order was small when the angular velocity was low and the sampling rate was large. However, as the angular velocity increased and the sampling rate decreased, the effect of integration order increased, i.e., obviously, the third-order KF resulted in better estimations than the first-order KF. When comparing the case where both transient matrix and process noise covariance matrix are applied to the corresponding order and the case where only the transient matrix is applied to the corresponding order but the process noise covariance matrix for the first-order is still used, both cases had almost equivalent estimation accuracy. However, in terms of the calculation cost, the latter case was more economical than the former, particularly for the third-order KF (i.e., the ratio of the former to the latter is 1.22 to 1).


Introduction
Applications of inertial measurement units (IMUs) consisting of accelerometers and gyroscopes have been exponentially increasing as these units can function as wearable motion sensors because of their sourceless property, i.e., they do not require external location-fixed sources [1][2][3][4][5][6]. Accurate attitude estimation based on an IMU is an important research theme, and therefore, several attitude estimation algorithms have been introduced in literature [7][8][9][10][11][12][13][14][15][16]. In spite of the varieties of the algorithms, the basic concept of estimation is common; that is, gyroscope signals are integrated to predict the attitude, and then accelerometer signals are used to prevent the drift error caused by the error accumulation associated with the integration.
The first prediction step is called strapdown integration as the sensors are rigidly strapped to a body that we want to track. With regard to the strapdown integration of the gyroscope signals, it should be noted that (i) a gyroscope actually provides discrete samples of the angular velocity, rather than continuous signals and (ii) the gyroscope signal is the summation of angular velocity and noise.
Related to (i), an integration scheme must be used to integrate the sampled signal and the choice of scheme is application-dependent. For short-timespan and low-accuracy applications, a low-order scheme may be sufficient. For more demanding applications, a higher-order scheme may be more appropriate, as discussed in [17]. It is simply natural that the higher-order scheme requires a higher computational cost. Once the integration order is chosen, the remainder of the Taylor series expansion is excluded from the integration equation. This approximation affects estimation accuracy to some extent. Related to (ii), when the gyroscope signals are integrated by an attitude estimation algorithm, errors in the signals propagate through to the calculated attitude. In the case of attitude algorithms using Kalman filters (KFs) which are the most popular cases, the noise-related terms form the process noise covariance matrices.
In literature, high-frequency angular motions-which are of interest in this study-are often discussed with "coning" or "noncommutivity rate" motions [18,19]. Coning motion (in which one or more axes of the system sweeps out a cone in space) is one particular input used to evaluate strapdown inertial navigation systems (INS) and attitude estimation algorithms. As coning is a demanding type of motion, it is admitted that algorithms operating satisfactorily in this environment will satisfy most other requirements [20][21][22][23][24][25]. However, most of these algorithms are in the context of aerospace fields (e.g., for spacecraft attitude determination) and are less related to biomedical or industrial applications that employ low-cost wearable sensors. Therefore, previous research has not investigated the effect of the difference in integration schemes in practical operating environments. The purpose of this research is to provide insights into the effect of the strapdown integration order and sampling rate on the attitude estimation accuracy for low-cost IMU applications.
In this paper, we compare and analyze IMU-based attitude estimation performance according to the strapdown integration order. Based on the attitude estimation KF introduced in [14], first-order, second-order, and third-order KFs are formulated with the first-order, second-order, and third-order strapdown integrations, respectively. Estimation performances of the three different KFs are evaluated under various conditions in terms of angular velocities and sampling rates as the transient matrix of KFs is the function of these two. In addition, we investigate the effectiveness of the exact implementation of the process noise covariance matrix of a KF by comparing the exact implementation to the approximated implementation. Finally, we discuss the estimation accuracy for each case, considering the calculation cost.

Methods
Measurements from the accelerometer (A) and gyroscope (G) sensors are respectively modeled as follows: where S g is the gravity vector, S a is the external acceleration, S ω is angular velocity, n A and n G are the measurement noises of the accelerometer and the gyroscope, respectively, which are assumed to be independent white Gaussian noise with zero mean and covariance matrices, respectively, are Σ A = σ 2 A I and Σ G = σ 2 G I. The left superscript S indicates that the vectors are observed in the sensor frame coordinate. In (1), the external acceleration S a is modeled as a first-order Markov chain stochastic process as in [14], i.e., where c a and c b are constant parameters and a w t−1 is the corresponding white Gaussian noise with zero mean. Both parameters c a and c b were set to 0.1 in this study. The attitude estimation KF proposed in [14] predicts the attitude in the time update (prediction step) using the angular velocity measured from the gyroscope and performs the measurement update (correction step) using the gravitational acceleration measured from the accelerometer. This paper is about the prediction step in which the selection of strapdown integration matters.

Strapdown Integration
The discrete-time model for the strapdown integration is [12]  where R is the short form of I S R, which is the rotation matrix of the sensor frame S with respect to the inertial frame I; ∆t is the sampling time; [ S ω×] represents the vector cross product of S ω = [ ω x ω y ω z ] T , i.e., Henceforth, [ S ω t−1 ×]∆t in (4) is simply denoted as A t−1 for convenience. By the Taylor expansion of the exponential term, (4) can be rewritten as Then, by transposing both sides of (6) and applying [a×] T = −[a×], (6) can be rewritten as Because the rotation matrix R contains the three unit column vectors of the inertial coordinate system expressed in the sensor coordinate system, i.e., T , the specific form of (7) for S Z is Again, this paper deals with the effect of the strapdown integration order on attitude estimation performance. Accordingly, depending on the strapdown integration order selected, (8) can be written as follows: where S Z 1,t , S Z 2,t , and S Z 3,t are the first-order, second-order, and third-order approximations of (8), respectively.

Transient Matrix and Process Noise Covariance Matrix
A linear KF can be defined by the following process and measurement models [14]: where x t is the state vector, Φ t−1 is the state transition matrix, w t−1 is the Gaussian process noise, z t is the measurement vector, H is the observation matrix, and v t is the Gaussian measurement noise. Covariance matrices of w t−1 and v t are Q t−1 and M t , respectively. The state vector is defined Note that S g is g × S Z where g is the norm of gravitational acceleration and attitude vector S Z is the sufficient information to calculate the attitude (i.e., roll and pitch) [12]. Therefore, the KF process model is where g Φ t−1 and g K t−1 vary according to strapdown integration orders. By substituting (2) into (9a-c) and applying [n G ×][n G ×] = 0, the process models in terms of y G,t and n G for the first-order, second-order, and third-order strapdown integrations are as follows, respectively: From (12), the g Φ t−1 's for the first-, second-, and third-order KFs are The process noise covariance matrices for the first-, second-, and third-order KFs are, respectively, where Q t−1 is defined as E w t−1 w T t−1 and E is the expectation operator. The attitude KF dealt in this paper has the measurement model based on the accelerometer measurement of (1) as follows, regardless of the integration order: where the measurement vector z t is y A,t , the observation matrix H is I I , the measurement noise v t is n A , and the measurement noise covariance matrix M t is Σ A [14].

Experimental Setup
For comparing the performance of KFs with respect to the strapdown integration order, an MTw IMU (from Xsens Technologies B.V., Enschede, The Netherlands) was used to provide input data for the aforementioned KFs. In addition, to investigate the estimation accuracy, an OptiTrack Flex13 3D optical tracking system (from NaturalPoint, Inc., Corvallis, OR, USA) was used to obtain the truth reference of the attitude based on the spatial positions of three markers on a plane (see Figure 1).
Because the angular velocity is one of the critical factors causing differences according to the integration order, three different experiments according to angular velocities were carried out. Each experiment repeated 30 times for Monte Carlo analysis. The means ± standard deviations of the averaged G y for each experiment are as follows: • The IMU measurement signals were delivered by the MTw at a sampling rate of 120 Hz. Then, to investigate the effects of the sampling rate in each test and for each integration order, the 120-Hz data were downsampled to 80-, 40-, and 20-Hz data by interpolation. Finally, a Monte Carlo analysis consisting of the 30 runs was performed. Because the angular velocity is one of the critical factors causing differences according to the integration order, three different experiments according to angular velocities were carried out. Each experiment repeated 30 times for Monte Carlo analysis. The means ± standard deviations of the averaged y G for each experiment are as follows: • Test A (slow): the average y G of nearly 1.47 ± 0.6 rad/s with the average maximum of 6.08 rad/s; • Test B (fast): the average y G of nearly 3.90 ± 0.6 rad/s with the average maximum of 13.70 rad/s; • Test C (vary fast): the average y G of nearly 6.18 ± 0.8 rad/s with the average maximum of 17.16 rad/s.
The IMU measurement signals were delivered by the MTw at a sampling rate of 120 Hz. Then, to investigate the effects of the sampling rate in each test and for each integration order, the 120-Hz data were downsampled to 80-, 40-, and 20-Hz data by interpolation. Finally, a Monte Carlo analysis consisting of the 30 runs was performed.

Case 1: Variation of Transient Matrix and Process Noise Covariance Matrix
Case 1 investigates the effect of the integration order on the estimation accuracy when both the transient matrix shown in (15) and the process noise covariance matrix shown in (16) are varied according to the selected order. In other words, the first-order KF uses g Φ 1 and Q 1 , the second-order KF uses g Φ 2 and Q 2 , and the third-order KF uses g Φ 3 and Q 3 . Table 1 shows the estimation results of roll, pitch, and attitude, in terms of the averaged root-mean-square error (RMSE) from 30 Monte Carlo runs, for different integration orders and for different sampling rates. The attitude error in Table 1 is the angle between the truth reference attitude vector from the optical tracking system ( S Z opt ) and the estimated attitude vector ( S Z est ). In Test A, although the RMSE decreases as the integration order increases, the integration order does not affect the accuracy as much as it does in Test B or C. As the sampling rate decreases, the RMSE increases slightly. In all cases, the RMSE was less than 2.2 • for roll and pitch, showing high estimation accuracy.
In Test B in which the angular velocity was higher than in Test A, RMSEs increased in all three KFs compared to those in Test A. For example, the attitude estimation RMSEs of first-, second-, and third-order KFs at 120 Hz were 3.11 • , 2.79 • , and 2.76 • , respectively. More importantly, the estimation performance in Test B was more affected by the selected integration order than that in Test A. The effect of the order on the RMSE was clearer in the lower sampling rate, i.e., the attitude estimation differences between the first-order KF and the third-order KF were 0.35 • at 120 Hz and 5.49 • at 20 Hz.
In Test C, in which the angular velocity was the highest among the three tests, the RMSEs also increased in all three KFs compared to those in Test B. For example, the attitude estimation RMSEs of first-, second-, and third-order KFs at 120 Hz were 5.26 • , 4.26 • , and 2.17 • , respectively. Note that IMU-based attitude determination under dynamic conditions is a type of "underdetermined" problems because the accelerometer signal used for the correction has two unknowns: the attitude and external acceleration. Therefore, particularly for highly dynamic test conditions like Test C, estimation performance can be seriously degraded even with a high sampling rate such as 120 Hz.
As shown in Figure 2 (from one trial out of 30 runs), in the case of a 40 Hz sampling rate for Test C, estimation errors from the first-order KF were significantly increased as time elapsed, compared to the second-order and third-order KFs. For example, the pitch estimation RMSEs from the first-order, the second-order, and the third-order KFs were 6.92 • , 4.22 • , and 3.19 • , respectively. However, when the sensor returned to the static condition, the estimation accuracy was recovered quickly in all KFs. The difference caused by the integration order increased, as the sampling rate decreased. For example, differences of the pitch estimation RMSEs from the first-order and third-order KFs were 0.75 • at 120 Hz, 1.47 • at 80 Hz, 3.79 • at 40 Hz, and 6.51 • at 20 Hz. This tendency can be observed in Test B as well.

KFs
Calculation Cost First-order KF 1 Second-order KF 1.12 Second-order′′ KF 1.04 Third-order KF 1.29 Third-order′′ KF 1.06 Therefore, when we need to choose a higher-order KF to deal with a high angular velocity and a small sampling rate and thus to maintain the estimation accuracy, it is more advantageous to change the transient matrix and to keep the process noise covariance matrix of the first-order integration. As seen in the three test results, the difference of the attitude estimation performances between the different order KFs became larger as A (= [ S ω×]∆t) was larger (i.e., as angular velocity S ω was larger and the time stepsize ∆t was larger and thus the sampling rate was smaller). Note that while all the KFs have the same correction procedure, more frequent correction procedures are applied to the cases of higher sampling rate than those of lower sampling rate. This means that the estimation differences shown in Table 1 do not come from the prediction procedure only.

Case 2: Variation of Transient Matrix Only
Case 2 investigates the effect of the integration order on the estimation accuracy when only the transient matrix shown in (15) is varied according to the selected order. In other words, the first-order KF uses g Φ 1 and Q 1 , the second-order KF uses g Φ 2 and Q 1 (which is referred to as second-order"), and the third-order KF uses g Φ 3 and Q 1 (which is referred to as third-order"). Table 2 shows the averaged RMSEs from the 30 Monte Carlo runs, for different integration orders and for different sampling rates. Similar to the results of Case 1, the difference of the attitude estimation performances between the different order KFs in Case 2 also became larger as the angular velocity was larger and the sampling rate was smaller. Note that in all the three tests, results from the second-order" and third-order" in Case 2 were almost the same as those from the second-order and third-order in Case 1, respectively. This shows that the effect of the process noise covariance matrix according to the integration order on the attitude estimation accuracy is negligible in most cases. It can be observed that some estimation errors from Case 1 were bigger than those from Case 2 (e.g., 9.44 • of attitude error from Case 1 versus 9.10 • from Case 2, in case of the third-order KF at 20 Hz of Test C). Such results occurred when low sampling rates (i.e., 20 Hz) and fast test conditions (i.e., Tests B and C) were applied and thus estimation accuracies were highly deteriorated.
In terms of computational costs for each approach compared to the first-order KF, the second-order" and third-order" in Case 2 were much more cost-effective than the second-order and third-order in Case 1, respectively, as the latter require the computation of the transient matrix as well as the process noise covariance matrix (see Table 3). In particular, the calculation time of the third-order" is 1.06 times higher than that of the first-order KF. On the other hand, the third-order KF requires 1.29 times more computation time than the first-order KF. However, when we compare the third-order" to the third-order, the former is 1.22 times cost-effective than the latter, whereas the former and the latter have almost the same estimation accuracy. Therefore, when we need to choose a higher-order KF to deal with a high angular velocity and a small sampling rate and thus to maintain the estimation accuracy, it is more advantageous to change the transient matrix and to keep the process noise covariance matrix of the first-order integration.

Conclusions
The attitude of an IMU relative to the global reference frame can be tracked by integrating the angular velocity signal obtained from the gyroscope in the IMU. In many low-cost and low-end IMU applications, a first-order integration scheme has been selected. With the improvement of the computer processing ability, the restriction on the calculation amount is reduced and the demand for the estimation accuracy is becoming stronger. Therefore, a higher-order integration scheme may be more easily chosen for low-cost applications than before.
In this paper, we investigated the effect of strapdown integration order and sampling rate on estimation accuracy under different test conditions. Experimental results showed that the effect of integration order was small when the angular velocity was low and the sampling rate was large. However, as the angular velocity increased and the sampling rate decreased, the effect of integration order increased, i.e., obviously, the third-order KF produced better estimations than the first-order KF. When comparing Case 1 (where both the transient matrix and the process noise covariance matrix are applied for the corresponding order) and Case 2 (where only the transient matrix is applied for the corresponding order but the process noise covariance matrix for the first-order is still used), both cases had almost equivalent estimation accuracy. However, in terms of the calculation cost, Case 2 is cheaper or more cost-effective than Case 1, particularly for the third-order (i.e., ratio of Case 1 to Case 2 is 1.22 to 1). Therefore, it can be concluded that Case 2 is superior to Case 1, overall.