State Derivation of a 12-Axis Gyroscope-Free Inertial Measurement Unit

Abstract

The derivation of linear acceleration, angular acceleration, and angular velocity states from a 12-axis gyroscope-free inertial measurement unit that utilizes four 3-axis accelerometer measurements at four distinct locations is reported. Particularly, a new algorithm which derives the angular velocity from its quadratic form and derivative form based on the context-based interacting multiple model is demonstrated. The performance of the system was evaluated under arbitrary 3-dimensional motion.

1. Introduction

Inertial sensors have been widely used in various applications, including motion detection [1], body state estimation [2–4], navigation [5–7], microsurgery [8], rehabilitation [9], etc. Traditionally a standard inertial measurement unit (IMU) comprised of 3-axis linear acceleration measurement by accelerometers installed at center of mass (COM) and 3-axis angular velocity measurement by rate gyros readily provides complete six degree-of-freedom (DOF) motion-related measurements spanning the 3-dimensional space. For highly dynamic systems which favorably have angular acceleration measurements, to the best of our knowledge there is no off-the-shelf product available. One widely adopted approach to derive this state is by differentiation of rate gyro signals, together with the filter technique. The other approach is based on the principle of Newtonian Mechanics, which relates linear acceleration, angular acceleration, and angular velocity in a memoryless manner. Because of this characteristic, derivation of the angular acceleration by only the inertial sensors seems to be a feasible method [10].

The gyroscope-free inertial measurement unit (GF-IMU) [11–13] is one of the more popular IMU methods to achieve this goal. Compared to the traditional IMU, the GF-IMU utilizing only accelerometers includes several features such as low-cost, easy calibration, being less affected by temperature variations, and a simple mechatronic setup. In general, the GF-IMU is capable of deriving linear acceleration, angular acceleration, and angular velocity. Because the latter two states have integrative/derivative relation, a GF-IMU comprised of 6-axis measurements is theoretically capable of yielding all three states (i.e., 9 scalar unknowns). One of the typical configurations of sensors is to have a 3-axis acceleration measurement at the COM and three 1-axis measurements on the principal axes. However, the iterative computation between the derived angular acceleration and the integrated angular velocity can possibly deteriorate the fidelity of these two states. Padgaonkar et al. proposed a 9-axis acceleration measurement system capable of deriving bounded linear and angular acceleration [14]. Chen et al. proposed a novel 6-axis system which yielded a bounded angular acceleration [15]. The system was carefully evaluated [16] and improved by adding a 3-axis acceleration measurement [17]. In general, due to the quadratic formulation of angular velocity in the rigid body dynamics, the derivation of this state in the 9-axis IMU faces the sign ambiguity problem [18]. This dilemma can be solved by comparing it to the estimated angular velocity which is integrated from the angular acceleration measurement [19] or by adding the redundant measurements to the IMU, for example, to increase the measurements to 12-axis [20]. Parsa et al. later developed an original all-accelerometer IMU which requires twelve 1-axis accelerometers mounted at specific locations on the surfaces of a cube. The system is capable of deriving all three states in which the angular velocity was derived through an optimization procedure from six measured inputs in the quadratic form [21]. Schopp et al. reported another novel 12-axis IMU which was constructed by twelve 1-axis accelerometers in different configurations and utilized an Unscented Kalman Filter (UKF) to yield all three states simultaneously [22].

Previously we had installed a 12-axis IMU composed of four 3-axis accelerometers at four distinct locations on the robot RHex [23], together with some custom-made leg sensors [24], to perform sensor data fusion for full body state estimation in this hexapod robot with dynamical gait [25]. Based on the rigid-body dynamics and matrix theory, the developed 12-axis IMU is theoretically capable of deriving all three states. However, limited available space on the RHex for sensor installation at that time constrained the configuration of the IMU far from the optimum level. Only the linear and angular accelerations were available for further analysis and no angular velocity developments were performed.

Here, we report on the state derivation and performance evaluation of the 12-axis IMU with optimal configuration in the sense of matrix operation, allowing the system to yield all three states. Particularly, a new algorithm which derives the angular velocity is reported. Basically, the state is estimated by the mixed signals from its quadratic form and derivative form based on the context-based interacting multiple models (IMM) [26]. The algorithm requires low computation power suitable for real-time derivation of the state. The proposed 12-axis IMU in its new configuration was tested under 3-dimensional random motion with various magnitudes, and its performance was evaluated by comparing to the results from the traditional IMU installed at the COM.

Section 2 briefly reviews the construction of the 12-axis IMU based on the analysis of rigid body dynamics. Section 3 describes the derivation of the angular velocity by the context-based IMM in detail. Section 4 reports the results of experimental evaluation, and Section 5 concludes the work.

2. Construction of the 12-Axis IMU

A brief review regarding construction of the 12-axis IMU is described in this section [25]. As shown in Figure 1, the acceleration vector, a_p, of a point, p, rigidly attached to an accelerating body frame B with origin o, in the inertial frame, W, is a function of the body’s angular velocity, ω, angular acceleration, ω̇, and translational acceleration of the body origin, a_o, represented by:

$${\mathbf{a}}_p={\mathbf{a}}_o+\dot{\mathbf{\omega }}\times {\mathbf{r}}_{\mathit{op}}+\mathbf{\omega }\times (\mathbf{\omega }\times {\mathbf{r}}_{\mathit{op}})$$

(1)

where r_op, the fixed position vector of point p relative to o, is assumed to be known. In general, the three body states (i.e., 9 scalar values) on the right hand side of Equation (1) are unknowns, including the COM translational acceleration, a_COM (usually equal to the origin of body frame, a_o), the body angular acceleration, and the angular velocity:

$$\begin{array}{l}{\mathbf{a}}_o={\mathbf{a}}_{\mathit{COM}}={\left[\begin{array}{ccc}{a}_x& a_y& a_z\end{array}\right]}^T\hfill \\ \dot{\mathbf{\omega }}={\left[\begin{array}{ccc}{\dot{\omega }}_x& {\dot{\omega }}_y& {\dot{\omega }}_z\end{array}\right]}^T\hfill \\ \mathbf{\omega }={\left[\begin{array}{ccc}{\omega }_x& {\omega }_y& {\omega }_z\end{array}\right]}^T\hfill \end{array}$$

(2)

With the quadratic representation of the angular velocity:

$${\mathbf{\omega }}_6(\mathbf{\omega })={\left[\begin{array}{llllll}{\omega }_x^2+{\omega }_y^2\hfill & {\omega }_x^2+{\omega }_z^2\hfill & {\omega }_y^2+{\omega }_z^2\hfill & {\omega }_x{\omega }_y\hfill & {\omega }_x{\omega }_z\hfill & {\omega }_y{\omega }_z\hfill \end{array}\right]}^T$$

(3)

Equation (1) appears to be linear with these 12 scalar unknowns:

$$\begin{array}{ll}{\mathbf{x}}_{\text{var}}\hfill & ={\left[\begin{array}{ccc}{\mathbf{a}}_o^T& {\dot{\mathbf{\omega }}}^T& {\mathbf{\omega }}_6{(\mathbf{\omega })}^T\end{array}\right]}^T\hfill \\ & ={\left[\begin{array}{llllllllllll}{a}_x\hfill & a_y\hfill & a_z\hfill & {\dot{\omega }}_x\hfill & {\dot{\omega }}_y\hfill & {\dot{\omega }}_z\hfill & {\omega }_x^2+{\omega }_y^2\hfill & {\omega }_x^2+{\omega }_z^2\hfill & {\omega }_y^2+{\omega }_z^2\hfill & {\omega }_x{\omega }_y\hfill & {\omega }_x{\omega }_z\hfill & {\omega }_y{\omega }_z\hfill \end{array}\right]}^T\hfill \end{array}$$

(4)

Presumably four 3-axis accelerometers are installed at point p_{j, j=1,2,3,4} with known r_{opj, j=1,2,3,4} :

$${\mathbf{r}}_{\mathbf{m}}={\left[\begin{array}{cccc}{\mathbf{r}}_{\mathit{op}1}^T& {\mathbf{r}}_{\mathit{op}2}^T& {\mathbf{r}}_{\mathit{op}3}^T& {\mathbf{r}}_{\mathit{op}4}^T\end{array}\right]}^T\text{with}\hspace{0.17em}{\mathbf{r}}_{\mathit{opj}}={\left[\begin{array}{ccc}{r}_{\mathit{opjx}}& r_{\mathit{opjy}}& r_{\mathit{opjz}}\end{array}\right]}^T j=1,2,3,4$$

and these accelerometers are oriented to measure accelerations in the directions along with three principal axes of the body coordinate, a_{pj, j=1∼4}:

$${\mathbf{a}}_{\mathbf{m}}={\left[\begin{array}{cccc}{\mathbf{a}}_{p1}^T& {\mathbf{a}}_{p2}^T& {\mathbf{a}}_{p3}^T& {\mathbf{a}}_{p4}^T\end{array}\right]}^T\text{with}\hspace{0.17em}{\mathbf{a}}_{\mathit{pj}}={\left[\begin{array}{ccc}{a}_{\mathit{pjx}}& a_{\mathit{pjy}}& a_{\mathit{pjz}}\end{array}\right]}^T j=1,2,3,4$$

(5)

a linear system with twelve equations and twelve unknowns is formed:

$$\begin{array}{ll}{\mathbf{a}}_{\mathbf{m}}\hfill & =\mathbf{S}\left({\mathbf{r}}_{\mathbf{m}}\right){\mathbf{x}}_{\text{var}}\hfill \\ & =\left[\begin{array}{cccccccccccc}1& 0& 0& 0& r_{1z}& -r_{1y}& 0& 0& -r_{1x}& r_{1y}& r_{1z}& 0\\ 0& 1& 0& -r_{1z}& 0& r_{1x}& 0& -r_{1y}& 0& r_{1x}& 0& r_{1z}\\ 0& 0& 1& r_{1y}& -r_{1x}& 0& -r_{1z}& 0& 0& 0& r_{1x}& r_{1y}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\end{array}\right]\hfill \end{array}{\mathbf{x}}_{\text{var}}$$

(6)

where S(r_m) is the 12 × 12 matrix and hereafter referred to as the “structure matrix”. The S(r_m) is the combination of four copies of Equation (1) with the dimensions 3 × 12. Due to the similarity of motion along with three principal axes, the structure of the 3 × 12 matrix is symmetric at a certain level. The first 3 × 3 matrix from the left side of S(r_m) is just an identity matrix and the second 3 × 3 matrix from the left side is the skew-symmetric matrix because of the cross product operator. The 3 × 6 matrix from the right side of S(r_m) is generated by the double cross product of the angular velocity term.

The unknown body states can now be derived by the matrix operation:

$${\mathbf{x}}_{\text{var}}=\mathbf{S}{\left({\mathbf{r}}_{\mathbf{m}}\right)}^{-1}{\mathbf{a}}_{\mathbf{m}}$$

(7)

Equation (7) reveals that the extraction of the desired state, x_var, now depends on the rank and numerical condition of the “structure matrix”, S(r_m), which is solely a function of the positions of accelerometers, r_m. Previously the numerical exploration pointed out that allocation of the four sensors shown in Figure 2(a) yields the best condition number of S(r_m), square root of 2. It indicates that this configuration is the most appropriate for matrix inversion [25], and the computation error induced by the matrix inversion is small. Therefore the new experimental setup shown in Figure 2(b) was built according to this configuration for the following analysis. Please note that the construction and inversion of the structure matrix S(r_m) only needs to be done once and can be computed offline after the positions and orientations of the accelerometers are determined.

3. Derivation of the Angular Velocity from the 12-Axis GF-IMU

Section 2 shows that in the real-time process the unknown body state x_var in the proposed 12-axis GF-IMU can be derived from the 12 sensed scalar accelerations multiplied with the inversed structure matrix shown in Equation (7). In x_var, the linear and angular accelerations are readily derived though the angular velocity in its exact form is still unsolved and requires further computation. In the current formulation, two sensed sources are available for this computation. One is from the angular acceleration, the 4th–6th terms of x_var shown in Equation (4), which is the derivative of the desired state. The other is from the last six terms of x_var, which is the quadratic form of the desired one. Because in the empirical setup the developed algorithm is executed by the commercial computers, the representation is in the discrete domain in the following sections.

3.1. The Primitive Derivation of the Angular Velocity

To derive the angular velocity from the available angular acceleration, ω_int, the trapezoid integration is the preferred method:

$${\mathbf{\omega }}_{\text{int}}: {\omega }_{i,j}={\omega }_{i,j-1}+\Delta t(\frac{{\dot{\omega }}_{i,j-1}+{\dot{\omega }}_{i,j}}{2}),i=x,y,z$$

(8)

where the subscripts j and j−1 represent consecutive two time stamps with time, Δt, in between. This method (hereafter referred to as the “integration method”) is simple and effective for a short duration, but is not suitable for a long period of time because the accumulated integration error quickly deteriorates the quality of the derived signal. Adding a bias term to Equation (8) may reduce the error; however, in general this compensation is not effective for motion not performed in specific patterns.

The angular velocity derived from its quadratic form (hereafter referred to as the “quadratic method”) has non-drift nature; however, the trade-off is the sign ambiguity problem, meaning to select the correct answer from multiple potential choices resulted from “square-root” computations. More specifically, assuming the first three quadratic terms shown in Equation (3) are chosen:

$$\begin{array}{ccc}{\omega }_x^2+{\omega }_y^2=a& {\omega }_x^2+{\omega }_z^2=b& {\omega }_y^2+{\omega }_z^2=c\end{array}$$

(9)

a solution derived from this method, ω_qua, can be represented as:

$$\begin{array}{ccc}{\omega }_{\mathit{qua},x}=\pm \sqrt{\frac{(+a+b-c)}{2}}& {\omega }_{\mathit{qua},y}=\pm \sqrt{\frac{(+a-b+c)}{2}}& {\omega }_{\mathit{qua},z}=\pm \sqrt{\frac{(-a+b+c)}{2}}\end{array}$$

(10)

which indicates that there are eight possible combinations because each scalar angular velocity has two possible solutions. When one or more values of the quadratic sums, a, b, and c, are very close to zero, the estimated angular velocity should be around zero, so the sign ambiguity problem vanishes and the number of combinations declines. However, in general situations when eight candidates appear, the selection process is required. Intuitively, the derived and readily-available angular acceleration shown in Equation (7) can be involved in this process. Without loss of generality, assuming ω_j−1 is the correctly estimated angular velocity at time stamp j − 1, the intuitive method to derive correct ω_j is to obtain the initial guess of ω_j by integration method shown in Equation (8), ω_int,j, and this guess is utilized as the reference to select one correct answer out of the eight candidates shown in Equation (10). More specifically, the L₂ norm can be utilized as the quantitative measure:

$${\Vert {\mathbf{\omega }}_{\mathit{qua}}-{\mathbf{\omega }}_{\mathit{int}}\Vert }_2$$

(11)

where the candidate with minimum value represents the correct choice. Figure 3(a) roughly sketches this computation process.

Practically, the quadratic method described above is likely to select an incorrect answer while the magnitude of the actual angular velocity, ω_act, is small. If the derived ω_j−1 at time stamp j − 1 is precise, the most likely cause of estimation error of ω_j at time stamp j “in the perfect world” is the process of trapezoid integration, which assumes the acceleration is constant during that time interval. As depicted in Figure 3(a), the actual motion pattern may vary in a very fast manner, and the angular acceleration derived from the 12-axis GF-IMU catches the instant value at the sampled time stamp because of the memoryless computation shown in Equation (7). The discrepancy between the instant and average accelerations during time stamp j results in the estimation error of the initial-guess, ω_int,j. This phenomenon in the traditional IMU (TIMU) which derives the angular acceleration by the differentiation of the angular velocity signal is even worse since the differentiation process introduces the noise and delay as shown in Figure 3(a). In the empirical world the situation is even more severe due to signal noises and accumulated digitization round-off errors during computation. For example, the estimated ω_j−1 at the j − 1 time stamp may already have certain estimation error, and the calculated angular acceleration and quadratic angular velocity shown in Equation (7) also have certain errors due to digitized matrix inversion and noisy sensor signals as depicted in Figure 3(b). Both empirical effects strongly affect the accuracy of 1-out-of-8 selection process, especially when the magnitude of the actual angular velocity ω_act is small, as depicted in Figure 3(c) where multiple choices of ω_qua may fall into the estimated ranges. In addition, because the estimation process is iterative, one incorrect estimate may badly affect the correctness of future estimates. Therefore it can be concluded that the quadratic method is suitable while the magnitude of quadratic sums, a, b, and c, are either very close to zero or large (i.e., not small).

In summary, neither one of the two methods is individually capable of yielding a correct estimation of the angular velocity. Because of the complementary characteristics between them, it is intuitive to fuse signals from these two methods to yield a better angular velocity estimate, ω_12-axis.

3.2. Context-Based Interacting Multiple Models

A better estimation of the angular velocity can be achieved by the adequate combination of the signals derived from the integration and quadratic methods. The process can be categorized within the domain of Interacting Multiple Models (IMM) [27,28], which generally calculates the accuracy of all models in a stochastic manner and mixes the estimated signals from all sources in a weighted manner to produce the correct estimate. Because executing the covariance of all models requires certain computation power as well as the performance of the models for specific scenarios may not be fairly judged by simple Gaussian assumptions, the context-based IMM [26] is adopted in the developed algorithm, which introduces the pre-selected contexts as the basic judgment for signal mixture from multiple models.

The development shown in the previous sub-section indicates that the quadratic method is effective while the magnitude of the angular velocity is either close to zero or very large. Therefore, two thresholds, T₁ and T₂, are selected as the contexts. T₁ is the boundary where the estimated angular velocity should be treated as zero, and T₂ is the boundary where the quadratic method is effective. These two contexts divide the range of quadratic sums, a, b, and c, into three sections as depicted in Figure 4. While 0 ≤ ω_i² + ω_j² ≤ T₁ as shown in Figure 4(a), referred to as the zero model, one can set the angular velocity to zero. While ω_i² + ω_j² ≥ T₂ as shown in Figure 4(c), referred to as the quadratic model, one can obtain the angular velocity by the quadratic method. While T₁ ≤ ω_i² + ω_j² ≤ T₂ as shown in Figure 4(b), referred to as the integration model, one adopts the integration method.

Each of the quadratic sums, a, b, or c, can reside in three possible sections shown in Figure 4, so there are twenty-seven possible combinations. Because the equations shown in Equation (9) are coupled, further categorization and treatment is detailed as follows.

Case 1. a, b, c ≥ T₂

When a, b, c ≥ T₂ as shown in Figure 5(a), ω_act is far away from zero along all three principal directions. In this case the estimated ω_12-axis is determined by the quadratic method:

$$\begin{array}{ccc}{\omega }_{12–\mathit{axis},x}={\omega }_{\mathit{qua},x}& {\omega }_{12–\mathit{axis},y}={\omega }_{\mathit{qua},y}& {\omega }_{12–\mathit{axis},z}={\omega }_{\mathit{qua},z}\end{array}$$

(12)

Case 2. T₁ ≤ a < T₂ and b, c ≥ T₂

When T₁ ≤ a < T₂ and b, c ≥ T₂ as shown in Figure 5(b), it is reasonable to conclude that ω_act,z is far away from zero and ω_act,x and ω_act,y are likely to have moderate magnitudes. Therefore both methods are utilized in this case:

$$\begin{array}{ccc}{\omega }_{12–\mathit{axis},x}={\omega }_{\mathit{int},x}& {\omega }_{12–\mathit{axis},y}={\omega }_{\mathit{int},y}& {\omega }_{12–\mathit{axis},z}={\omega }_{\mathit{qua},z}\end{array}$$

(13)

Similarly, both T₁ ≤ b < T₂, a, c ≥ T₂ and T₁ ≤ c < T₂, a, b ≥ T₂ are within this case.

Case 3. 0 ≤ a < T₁ and b, c ≥ T₂

When 0 ≤ a < T₁ and b, c ≥ T₂ as shown in Figure 5(c), it can be concluded that ω_act,z has a large magnitude but ω_act,x and ω_act,y are close to zero. Therefore only the former requires computation by the quadratic method:

$$\begin{array}{ccc}{\omega }_{12–\mathit{axis},x}=0& {\omega }_{12–\mathit{axis},y}=0& {\omega }_{12–\mathit{axis},z}={\omega }_{\mathit{qua},z}\end{array}$$

(14)

Similarly, both 0 ≤ b < T₁, a, c ≥ T₂ and 0 ≤ c < T₁, a, b ≥ T₂ are within this case.

Case 4. 0 ≤ a, b, c < T₁

When 0 ≤ a, b, c < T₁ as shown in Figure 5(d), it is reasonable to set all components to zero:

$$\begin{array}{ccc}{\omega }_{12–\mathit{axis},x}=0& {\omega }_{12–\mathit{axis},y}=0& {\omega }_{12–\mathit{axis},z}=0\end{array}$$

(15)

Besides the four cases shown above, there are nineteen combinations left undetermined. Because there is no clear trend to judge the adequateness of the quadratic method in these combinations, the estimated angular velocity ω_12-axis is derived from the integration method (i.e., ω_12-axis = ω_int). If the computed a, b, or c appears in an unreasonable less-than-zero value due to empirical computation error, the estimated angular velocity ω_12-axis is also derived from the integration method (i.e., ω_12-axis = ω_int, the quadratic sums are not utilized in the computation in this time stamp).

The proposed estimation shown in the previous paragraph has the “hard switching” nature. This implies an unreasonable situation because the trustworthiness of the models has a sharp boundary. Technically, suddenly switching the estimation from one method to another also introduces a signal discontinuity problem. Therefore the “soft switching” technique is adopted, which defines the probability of each model in a continuous manner as shown in Figure 6. When ω_i² + ω_j² is close to T₁, ω_12-axis is designed to be the linear combination of zero and ω_int :

$${\mathbf{\omega }}_{12–\mathit{axis}}=(1-\frac{{\omega }_i^2+{\omega }_j^2-T_1}{h_1\Delta })\mathbf{0}+(\frac{{\omega }_i^2+{\omega }_j^2-T_1}{h_1\Delta }){\mathbf{\omega }}_{\mathit{int}}$$

(16)

where Δ = T₂ − T₁. When ω_i² + ω_j² is close to T₂, ω_12-axis is designed to be the linear combination of ω_int and ω_qua :

$${\mathbf{\omega }}_{12–\mathit{axis}}=(\frac{T_2-({\omega }_i^2+{\omega }_j^2)}{h_2\Delta }){\mathbf{\omega }}_{\mathit{int}}+(\frac{({\omega }_i^2+{\omega }_j^2)-(T_2-h_2\Delta )}{h_2\Delta }){\mathbf{\omega }}_{\mathit{qua}}$$

(17)

The h₁ and h₂ are the percentages of the overall range to be utilized for linear combination of models (i.e., soft switching). Smaller values of h₁ and h₂ represent sharper switching, and larger values represent smoother transition. The h₁ and h₂ are set around 2% in the empirical evaluation.

3.3. Brief Discussion

The algorithm reported in the previous sub-section utilizes the first three terms shown in Equation (3) to recover the angular velocity. The last three terms can also be utilized to perform the quadratic inversion. Assuming the quadratic multiplication is labeled as:

$$\begin{array}{ccc}{\omega }_x{\omega }_y=d& {\omega }_x{\omega }_z=e& {\omega }_y{\omega }_z=f\end{array}$$

(18)

and the angular velocity can be derived as:

$$\begin{array}{ccc}{\omega }_x=\pm \sqrt{\frac{\mathit{de}}{f}}& {\omega }_y=\pm \sqrt{\frac{\mathit{df}}{e}}& {\omega }_z=\pm \sqrt{\frac{\mathit{ef}}{d}}\end{array}$$

(19)

Because of the multiplication characteristics there are only two candidates instead of eight. This reveals that if one sign of ω_{i, i=x,y,z} is selected, the other two will be determined. However, empirically the sign determination is tricky and no obvious model can be constructed. Thus this approach is not adopted in the development. Another method is to construct the complete squares by utilizing all six terms shown in Equation (3):

$$\begin{array}{cccccc}{\omega }_x^2+{\omega }_y^2=a& {\omega }_x^2+{\omega }_z^2=b& {\omega }_y^2+{\omega }_z^2=c& {\omega }_x{\omega }_y=d& {\omega }_x{\omega }_z=e& {\omega }_y{\omega }_z=f\end{array}$$

(20)

and angular velocity can be derived as:

$$\begin{array}{l}{\omega }_1+{\omega }_2=\pm \sqrt{{\left({\omega }_1+{\omega }_2\right)}^2}=\pm \sqrt{\left(a+2d\right)}\hfill \\ {\omega }_1+{\omega }_3=\pm \sqrt{{\left({\omega }_1+{\omega }_3\right)}^2}=\pm \sqrt{\left(b+2e\right)}\hfill \\ {\omega }_2+{\omega }_3=\pm \sqrt{{\left({\omega }_2+{\omega }_3\right)}^2}=\pm \sqrt{\left(c+2f\right)}\hfill \end{array}$$

(21)

Though the sensed signals shown in Equation (7) are utilized in a more thorough manner, the strategy of setting adequate contexts is also not clear in this case. For example, the advantage of the zero model does not exist in this case because the setting of ω_x + ω_y = 0 can only reveal that ω_x and ω_y are in opposite sign with no magnitude information. Thus this approach is not adopted in the development either.

4. Experimental Evaluation

The experimental apparatus shown in Figure 2(b) was utilized for experimental evaluation of the proposed 12-axis system. The required sensory measurements were obtained from four 3-axis accelerometers (ADXL330, ±3g, Analog Device) installed at the specific configuration shown in Figure 2(a). In addition, a traditional IMU composed by one 3-axis accelerometer (ADXL330) and three 1-axis rate gyros (ADXRS610, ±3000/s, Analog Device) was also mounted at the COM for performance comparison. A real-time embedded control system (sbRIO-9632, National Instruments) running at 500 Hz was in charge of sensor signal collection. All of the analog input channels of the sbRIO have ±10 V range and 16-bit A-to-D resolutions. Random motions with varied magnitudes were applied to the experimental apparatus during experiments and the following analysis was based on the measured data.

4.1. Selection of Contexts T₁ and T₂

The context T₁ represents the boundary which sets the estimated angular velocity ω_12-Axis at zero. It is not reasonable to set a large T₁ as it would force ω_12-Axis to be zero when it is not. On the other hand, a very small T₁ yields very little data that qualifies for this criterion. Empirically it is determined by the noise level of the sensors as well as the precision of the digitized computation. T₁ is set around 0.1 in the experiments.

The context T₂ determines the magnitude level where the angular velocity can be effectively determined by quadratic methods instead of the integration method. Therefore small T₂ easily yields the wrong selection from the eight candidates. Large T₂ forces the data to be computed by integration, and the data drift appears when the time duration of the angular velocity computed in this method is long. Therefore a study on how to correctly choose the right T₂ is performed and detailed as follows.

Figure 7 plots the typical Root-mean-squared Error (RMSE) vs. T₂ based on one of the experimental data, where the RMSE is the comparison between the estimated angular velocity, ω_12-Axis, and that measured by the traditional IMU, ω_TIMU. The RMSE shown in the plot is the summed result of its three scalar components. The plot indicates that the RMSE is relatively large when T₂ is small, when the quadratic method is over trusted. It also indicates that RMSE is relatively large when T₂ is large. In this setting most of the estimates were done by the integration method and the data drift was observed. The wide and flat bottom of the curve shown in Figure 7 is also observed in other data sets, which indicates that there exists a wide selectable range of T₂ values which yields similar performance, as the best RMSE happened at T_{2 min}. For example, if the acceptable RMSE is bounded by an extra 20% of the best RMSE, the selectable range of T₂ is spanned from 3 to 9.

Figure 8 plots the variation of T_{2 min} (blue circle) vs. the average level of the motion, which is quantitatively defined as the summation of the quadratic sums, a + b + c. Instead of defining the level of motion directly in the angular velocity, quadratic sums are utilized since these sums are available right after the multiplication of the inversed matrix to the sensory inputs shown in Equation (7). Because the errors resulted from the matrix inversion and noises due to empirical sensor readings are usually scaled with the magnitude of the signals, the selected T₂ should increase as the magnitude of motion increases. The blue linear trend line with positive slope also confirms this phenomenon. The plot also reveals that the tolerable 10% or 20% increase of RMSE intersects with the linear trend line. Because the lengths of 20% lines are large and the slope of the trend line is small, the computation error of a + b + c has very little effect on the quantitative measure of the trend line. Thus the adequate T₂ can be easily obtained according to the equation of the trend line when the quadratic sums, a, b, and c, are given. This plot suggests that the selection of T₂ can be achieved by the given quadratic sums and the trend line, and the RMSE comparison test which requires the gyroscope input shown in Figure 7 is not necessary. The selected T₂ is fixed for the followed real-time estimation.

In order to quantitatively evaluate the usability of the trend line, instead of using T_{2 min} as the context, Figure 9 plots the percentage error of the estimated angular velocity versus T_{2 trend}, which is the selected T₂ calculated from the trend line with given quadratic sums. Percentage error is calculated as the ratio of the RMSE to the maximum magnitude at that experiment trial, where the RMSE is the comparison between the estimated angular velocity, ω_12-Axis, and that measured by the traditional IMU, ω_TIMU. Figure 9 indicates that the computed T₂ from trend line performs adequately; the percentage errors are mostly around 12% or less.

4.2. Performance of the State Derivation from the 12-Axis GF-IMU

In the experimental evaluation the apparatus was moved arbitrarily in the 3-dimensional space; thus the linear and angular accelerations along all three principal axes could be induced in the test. Before the sensor readings were imported into Equation (5), the raw accelerometer readings were filtered by Chebyshev filters and gravity-compensated by the readings of the 2-axis inclinometer.

Table 1. The RMSE Between the traditional IMU and the 12-axis IMU while the apparatus moved arbitrarily in the 3-dimensional space.

Linear acceleration (m/s2)			Angular acceleration (rad/s2)			Angular velocity (rad/s)
ax	ay	az	ω̇x	ω̇y	ω̇z	ωx	ωy	ωz
0.1359	0.0933	0.1296	0.4985	0.7691	1.4188	0.3953	0.2301	0.2593

lists the statistical summary of the experiments where the RMSEs were the comparisons between the estimated state of ω_12-Axis and the measured state from the traditional IMU installed at the COM, ω_TIMU. The angular acceleration of the traditional IMU is obtained by differentiation of the angular velocity, followed by filtered with a Chebyshev filter. Figure 10 plots one typical result of the experiment. Figure 10(a) confirms that though the sensors of the proposed 12-axis GF-IMU (thin red solid lines) are not located at the COM, they can indeed recover the COM linear acceleration. Figure 10(b) shows that the angular acceleration can also be correctly derived by the 12-axis GF-IMU. Figure 10(c) reveals that though several unmatched sections exist between the 12-axis GF-IMU and the traditional IMU readings, the proposed algorithm in general can indeed recover the angular velocity along with all three principal directions. The discrepancy either resulted from (i) the accumulated integration error where the magnitude of the quadratic sums fell into the integration model for too long or (ii) the incorrect selection of the angular velocity in the quadratic model.

Figure 11 compares the performance of three different methods: integration method, quadratic method, and the proposed method. Figure 10(c) reveals that the angular velocity derived from the proposed method matches closely to the readings from the traditional IMU, so in Figure 11 the performance of the latter one is skipped for clear presentation. Figure 11 shows that the angular velocity derived from the integration method drifts over time as expected. In contrast, the angular velocity derived from the quadratic method is bounded. However, the 1-out-of-8 selection process of the quadratic method is likely to select an incorrect answer while the magnitude of the signal is small. In addition, because the estimation process is iterative, one incorrect estimate may badly affect the correctness of future estimates until at certain moment the correct selection moves the estimates back to the right track. Figure 11 indicates that the proposed method with right mixture between the integration and the quadratic methods yields the adequate estimation.

Figure 12 shows the timings where the switching between two methods takes place. The 15-sec data had 7,500 sampled data points, and the proposed method switched around 100 times.

5. Conclusions

We have investigated a 12-axis inertial measurement unit that utilizes four 3-axis linear acceleration measurements from accelerometers installed at four distinct locations. We have developed a new algorithm which derives the angular velocity by mixing the signals from its quadratic form and its derivative form via the context-based interacting multiple models. The performance of the system was evaluated while the system was under arbitrary 3-dimensional motion. By adequately-choosing two contexts, the angular velocity can indeed be recovered. In the meantime, the linear and angular accelerations are correctly estimated as well, which confirmed that the COM acceleration state can be derived even though the sensors are not installed at that specific spot.

We are in the process of investigating a sensor fusion scheme of the reported system with other position and orientation sensors with the intention of constructing an observable system capable of accurate full body state estimation for analysis of dynamic locomotion in legged robots.