Angular Rate Estimation Using a Distributed Set of Accelerometers

A distributed set of accelerometers based on the minimum number of 12 accelerometers allows for computation of the magnitude of angular rate without using the integration operation. However, it is not easy to extract the magnitude of angular rate in the presence of the accelerometer noises, and even worse, it is difficult to determine the direction of a rotation because the angular rate is present in its quadratic form within the inertial measurement system equations. In this paper, an extended Kalman filter scheme to correctly estimate both the direction and magnitude of the angular rate through fusion of the angular acceleration and quadratic form of the angular rate is proposed. We also provide observability analysis for the general distributed accelerometers-based inertial measurement unit, and show that the angular rate can be correctly estimated by general nonlinear state estimators such as an extended Kalman filter, except under certain extreme conditions.


Introduction
A distributed accelerometers-based inertial measurement unit (IMU) uses only accelerometers to compute a specific force and angular rate. It has several advantages over the conventional gyroscope-equipped IMU, including robustness, easy calibration, low cost, and less power consumption. Furthermore, it can also function as an angular acceleration sensor, and thus it has applications in highly dynamic systems area, as well as low-cost, medium performance inertial navigation systems.

OPEN ACCESS
Development efforts for a distributed accelerometers-based IMU have been going on for over 30 years [1]. The research has mainly focused on the design of an optimal accelerometer configuration capable of overcoming the inherent calculation error that increases with time [2]. Although it was known in theory that a minimum of six accelerometers are required for a complete description of a rigid body motion, such six accelerometer schemes were not realized until Chen [3] proposed a cube-type IMU, which has one accelerometer at the center of each surface of a cube and its sensing direction is along the respective surface diagonal. This system was carefully evaluated and it was proved by Tan [4] that any six accelerometer configuration can be converted to the cube-type IMU with same computational simplicity. However, a six accelerometer-based IMU requires extra integration to obtain the angular rate, and thus a distributed accelerometers based inertial navigation system will have angular rate estimates that rapidly diverge, where the divergence rate is an order of magnitude greater than that of a gyroscope-equipped system. Park [5] proposed a redundant accelerometer-aided IMU with nine accelerometers by adding an extra three-axis accelerometer to the cube-type IMU to remove the extra integration using an extended Kalman filter, but his scheme suffers from an observability problem if one of the angular rate terms is zero.
The minimum number of accelerometers necessary to directly calculate the angular rate for a rigid body motion is 12 [6]. The accelerometer-based IMU with a minimum number of 12 accelerometers allows one to compute the angular rate without using the process of integration of the output of the accelerometers, and thus it can significantly improve the performance of the IMU. This scheme was studied further by Parsa [7] and Lin [8], and experimentally validated by Cappa [9]. However, the main problem with the 12 accelerometer-based IMU is that the angular rate is present in the system equations in its quadratic form, and therefore, it is difficult to extract the magnitude of the angular rate in the presence of accelerometer noises, and even worse, it is not possible to determine the direction of the angular rotation using only the quadratic form of the angular rate. In the literature, various methods have been proposed to solve this problem. Parsa [7] used a nonlinear least-square optimization procedure to estimate the angular rate from six measured quadratic forms of the angular rate. Cardou [10] reviewed various schemes available in the literature, and proposed a new algorithm for the estimation of the angular rate from the centripetal components of the acceleration measurements, while the sign of the angular rate was simply chosen by comparing the estimate with that of the integration of angular acceleration. Their schemes focused on extracting only the magnitude of angular rate in the presence of accelerometer noises. Continuous research efforts have been carried out to fuse the angular acceleration and the quadratic form of the angular rate. Schopp [11] applied an unscented Kalman filter to the entire 12 system equations to determine all the kinematic states such as specific force, angular rate and angular acceleration at the same time, but the observability was not proven. Lu [12] developed a new algorithm which derives the angular rate by mixing the signals from its quadratic form and its derivative form via context-based interacting multiple models, but his scheme only utilized the first three of the six quadratic angular rate terms, and the adequate context was set heuristically.
In contrast to those solutions, we propose in this paper an extended Kalman filter scheme to aid the integration of angular acceleration by six quadratic terms of angular rate in order to correctly estimate both the direction and magnitude of the angular rate. Compared to the previous works, we apply a Kalman filter to estimate angular rate only, since the specific force and angular acceleration can be computed algebraically. For an extended Kalman filter setup, the 4th-6th variables among the computed 12 kinematic variables in the system equation are designated as a state propagation equation, and the last six variables are designated as a measurement equation, which is the approach originally introduced in [5] by the authors. Furthermore, in this paper, we pay special attention on exploring the fundamental limit of the general distributed set of accelerometers based IMU using observability analysis, and formally show that the angular rate can be correctly estimated by general nonlinear state estimators such as the extended Kalman filter, except under certain extreme conditions.

Kinematics of the Distributed Accelerometers Based IMU
Consider the inertial frame {i} and a point k fixed in a rigid body moving in space, to which a body frame {b} is attached, as shown in Figure 1.
is the position vector from the center of the inertial frame to the center of the body frame, is the position vector from the center of the inertial frame to a point k, and is the position vector from the center of the body frame to the point k. Then, the acceleration of the point k with respect to the inertial frame is given by: where is the specific force at the point k and g i is the gravitational acceleration, and both are represented in the inertial frame {i}. Vector is represented by in the body frame {b}, and is a direction cosine matrix that takes frame {i} to frame {b}. The term is the angular rate of frame {b} with respect to the frame {i}, represented in frame {b}, and [ ×] is a cross-product matrix of the angular rate = [ω 1 ω 2 ω 3 ] T , which is given by: If an accelerometer is rigidly attached at point k with the sensing direction , the output a k of the accelerometer is given by: The output of the accelerometer is directly related with the specific force at the center of the body frame {b}, , the rigid body angular acceleration , whose components appear in the skew-symmetric elements of [ ×] as follows: (3) and the angular rate , whose components appear as quadratic products in the elements of [ ×] 2 as follows: Through simple algebraic manipulations, it is easily shown that Equation (2) can be expressed with 12 kinematic variables as follows: [ where in f i 's, r i 's, and s i 's (for i = 1, 2, 3) denote the components of , and , respectively. Now suppose that N accelerometers are distributed in the body frame, and if we collect the accelerometer outputs in a vector form, we have: where the N × 12 matrix J is called a configuration matrix, and it consists of the constant parameters about the positions and the sensing directions of the accelerometer array.
If the configuration matrix J has a left inverse matrix, then it is possible to calculate 12 kinematic variables as follows: where J + = (J T J) −1 J T is left inverse matrix of J, for which to exist, a minimum of 12 accelerometers is necessary. Equation (9) is called the system equation. For a configuration of 12 accelerometers, the left inverse matrix J + becomes the inverse matrix J −1 , and the main requirement for the placement of the accelerometers is that the configuration matrix is invertible. Several accelerometer array configurations have been proposed in the literature, and  Additional requirements can be imposed to determine the positions and the sensing directions of the accelerometers, since we have 60 parameters for an array of 12 accelerometers to be determined. For example, we can try to find an accelerometer array configuration yielding the minimum influence of the accelerometer noise on the variances of the kinematic variables, since it depends on the accelerometer configuration as follows: where X y and X A are the covariance matrices of the kinematic variables and the accelerometer array, respectively, and is the variance of the single accelerometer. The second term of Equation (10) is valid under the assumption that the accelerometer array is i.i.d. For another example, Lin [8] studied the accelerometer array configuration producing a minimum condition number of J, and pointed out that placement of the four 3-axis accelerometers shown in Figure 2(b) yields the best condition number of J, which indicates that this configuration is the most appropriate for matrix inversion, although the inversion of the configuration matrix J can be computed offline only once.

Estimation of the Angular Rate
The 12 kinematic variables can be computed from Equation (9), and the specific force at the center of the body frame, angular acceleration, and the magnitude of the angular rate can be obtained. Therefore, a distributed accelerometer-based IMU can perform the same function as a conventional gyro-equipped IMU, except for the direction of the angular rotation, whose magnitude is directly measured by gyroscopes in a conventional gyro-equipped IMU. It is not straightforward to determine the sign of the angular rate using a distributed accelerometer-based IMU since the angular rate is present in its quadratic form, as shown in Equation (6).
Several methods have been proposed to resolve this sign ambiguity. Some methods resort to additional sensors such as the conventional low-cost gyroscopes or GPS to determine the correct sign of the angular rate. However, the common approaches compare the estimates with those obtained from the integration of the angular acceleration. In this paper, we use an extended Kalman filter to fuse two sensed sources, the angular acceleration and the quadratic form of the angular rate, and thus to estimate the sign and magnitude of the angular rate at the same time.
For an extended Kalman filter setup, the angular accelerations, which are the 4th-6th terms of y in Equation (6), are designated to establish the angular rate propagation equation, and the quadratic forms of the angular rate, which are the last six terms of y, are designated to construct angular rate measurement equation.
Then, the angular rate propagation and measurement equations are written as:

Observability Analysis
Observability denotes the ability to uniquely determine the states of a dynamic system from its measurements and it is the indication for the feasibility of design of a state estimator such as a Kalman filter [13]. Observability can be examined by the rank of the corresponding observability matrix, and the observability matrix of nonlinear systems (11) and (12) is given as follows: where W 0 is a 18 × 3 matrix, and stands for a gradient operator with respect to , which is given by: As seen from Equation (13), the observability property can be changed by the angular acceleration as well as the angular rate itself to be measured. There are three cases to be considered for the observability analysis.
In the first case, wherein at least one of angular acceleration terms is not zero, the rank of observability matrix is three and thus the system is observable for every angular rate including the origin of angular rate vector space, i.e., = [0 0 0] T , which means that the system is globally observable in the entire angular rate vector space.
In the second case, in which all of angular rate terms are constant, the system becomes static but it is still observable if at least one of angular rate terms is not zero since the rank of observability matrix is three. However, the observability condition fails at the origin of angular rate vector space, which means that the system is locally observable in the angular rate vector space. Indeed, it is not possible to determine the direction of angular rate using only quadratic terms of the angular rate, since the quadratic terms do not give a unique solution. Instead, it yields two possible solutions, each of which having equal magnitude but opposite direction in the angular rate vector space, as shown in Figure 3. Thus, if an a priori estimate of angular rate is closer to the true angular rate than to the false angular rate, then the estimate converges to the true value. However, if the a priori estimate of angular rate is closer to the false angular rate, then the estimate goes to the false value. Thus, if there is no information on the angular rate at a previous time, it cannot be guaranteed for the angular rate to be correctly estimated. In the general application there are certain instants where the second case may occur. However, this is not a problem because the recursive predictor-corrector scheme of Kalman filter makes it possible to continuously incorporate an a priori estimate of angular rate, which must be close to the true angular rate, into the filtering process. Instead, the extremely severe situation is the case in which all of the angular rates are constant from the initial time when the system starts to work, because there is no available information for the initial guess on the direction of a rotation. Although this hardly occurs in real world applications, this situation may be avoided by flipping the sign of angular rate estimate where necessary when the divergence of the filter is detected.
The last case is a special case in which all of angular rate and angular acceleration terms are zero. Since all elements of observability matrix are zero, the system becomes completely unobservable. However, although the system is completely unobservable, this case can be treated easily since we can be informed that all angular accelerations and rates are zero from the designated measurement equation.

Extended Kalman Filter and Its Modification
The angular rate propagation Equation (11) and measurement Equation (12) are re-written here in discrete-time form: where the subscript k denotes the time update sequence and ΔT is the sampling interval. The terms w(t k ) and v(t k ) are included to account for noise in the computations of the angular acceleration and the quadratic forms of the angular rate with Equation (9). They originate from noises in the accelerometer outputs, which are assumed to be uncorrelated zero-mean white with the known power spectral densities. The covariance matrices Q and R of w(t k ) and v(t k ) can be computed using Equation (10).
An extended Kalman filter can be designed using these two equations. The time update of the angular rate estimate and its covariance matrix are given by: where (t k+1 |t k ) and P(t k+1 |t k ) are the a priori estimates of angular rate and error covariance. The time update portion of the extended Kalman filter gives a prediction of the angular rate at time t k+1 , along with the associated error covariance P. The extended Kalman filter gain, angular rate estimate and its covariance are updated as follows [14]: The measurement update provides a correction based on the measurement z at time t k+1 to yield the a posteriori estimate and its covariance. For the case in which all of both angular acceleration and quadratic angular rate terms lie under the certain thresholds and can be considered as zero, the measurement update portion can be modified as follows: The thresholds for angular acceleration and quadratic angular rate terms can be chosen by several methods. One method may be using the standard deviation of their measurement noises, and another method may be using a moving average of the measurements to obtain smoother transition.

Simulation Results
To evaluate the performance of the proposed scheme for the angular rate estimation, computer simulations are performed. We use the 12 accelerometer-based IMU shown in Figure 2(b). The IMU cube length l is 10 cm, and the position of the accelerometer array is determined as follows: and the sensing direction is determined as follows: The noise incorporated in the simulation is white, and the noise associated with each accelerometer measurement is assumed to have the same density of 100 μg/√Hz. The other errors such as scale factor error, bias, and the configuration (position and sensing direction) error of the accelerometer are assumed to be identified and compensated by some calibration scheme beforehand. The time update of the IMU algorithm is 100 Hz. Simulations are performed for three different cases as follows.

Case 1: Arbitrary Trajectory
The arbitrary trajectory shown in Figure 4 is used for the simulations. Figure 4 shows that the specific force and angular acceleration converge to their true values. Since the calculation of specific force and angular acceleration uses measurements obtained directly from the accelerometers, it never fails to identify them. In this simulation, the trajectory to be followed is sinusoidal. Thus, the condition of observability is satisfied and the extended Kalman filter is expected to estimate the angular rate correctly, as discussed previously. Figure 5 shows that the angular rate is correctly estimated by the extended Kalman filter.  Figure 6 shows the trajectory of angular rate ω 3 , used for the simulations. Other angular rate and angular acceleration terms are assumed to be zero, as shown in Figure 7. The extended Kalman filter correctly estimates the angular rate including the direction of it, which is expected by the observability analysis. Note that no drift of the estimates occurs when no angular rate is present, thanks to the modified measurement update portion.

Case 3: Constant Angular Rate
The extremely severe situation for distributed an accelerometer-based IMU is the case in which all of the angular rates are constant from the initial time when IMU is on, because there is no available information on the direction of rotation. According to the observability analysis, the extended Kalman filter cannot guarantee that the angular rate estimate will converge toward the correct values, unless the initial angular rate estimate is closer to the true angular rate than to the false angular rate.
The simulation was performed with two different initial estimates within the extended Kalman filter. Figure 8 shows the simulation results when the initial angular rate estimate is closer to the true angular rate. As discussed in the observability analysis, the estimate converges to the true angular rate. Figure 9, on the other hand, shows that the estimate converges to the false angular rates, which has the same magnitude but opposite direction and thus produces same values of quadratic angular rate terms. The two results agree well with our analytic arguments. The situation shown in Figure 9 can be avoided when the angular acceleration starts to vary as shown in Figure 10, where the trajectory of angular acceleration follows a sinusoidal pattern at 65 s. As shown in Figure 11, although initially the estimate converges to the false angular rate, eventually it converges to the true value. This is achieved by flipping the sign of angular rate estimate when the divergence of the filter is detected around at 67 s. In this simulation, the divergence of the filter is detected based on the magnitude of the measurement prediction covariance of the filter as shown in Figure 12. The measurement prediction covariance, S(t k ), is given by following equation [14]:

Conclusions
Although a typical 12 distributed accelerometer-based IMU allows for the computing of the magnitude of an angular rate without using the integration of the outputs of accelerometers, it is difficult to determine the direction of a rotation because the angular rate is present in its quadratic form within the IMU system equations. To tackle this inherent disadvantage, we proposed an extended Kalman filter scheme to aid the integration of angular acceleration by six quadratic terms of angular rate in order to correctly estimate both the direction and magnitude of the angular rate. We also provided observability analysis for the general 12 accelerometer-based IMU, and showed that the angular rate can be correctly estimated by general nonlinear state estimators such as an extended Kalman filter, except under certain extreme conditions. The performance of the proposed scheme for the angular rate estimation was evaluated using computer simulations under three special cases. The simulation results agree well with our analyses, and show that the proposed scheme correctly estimates specific force, angular acceleration, and the direction and magnitude of the angular rate as well.