A Multi-Model Combined Filter with Dual Uncertainties for Data Fusion of MEMS Gyro Array

The gyro array is a useful technique in improving the accuracy of a micro-electro-mechanical system (MEMS) gyroscope, but the traditional estimate algorithm that plays an important role in this technique has two problems restricting its performance: The limitation of the stochastic assumption and the influence of the dynamic condition. To resolve these problems, a multi-model combined filter with dual uncertainties is proposed to integrate the outputs from numerous gyroscopes. First, to avoid the limitations of the stochastic and set-membership approaches and to better utilize the potentials of both concepts, a dual-noise acceleration model was proposed to describe the angular rate. On this basis, a dual uncertainties model of gyro array was established. Then the multiple model theory was used to improve dynamic performance, and a multi-model combined filter with dual uncertainties was designed. This algorithm could simultaneously deal with stochastic uncertainties and set-membership uncertainties by calculating the Minkowski sum of multiple ellipsoidal sets. The experimental results proved the effectiveness of the proposed filter in improving gyroscope accuracy and adaptability to different kinds of uncertainties and different dynamic characteristics. Most of all, the method gave the boundary surrounding the true value, which is of great significance in attitude control and guidance applications.


Introduction
The advent of micro-electro-mechanical system (MEMS) gyroscopes have added motion sensing in consumer devices and special military applications [1][2][3]. This is because it has many excellent properties, including being lightweight and having a compact size, low power consumption, low cost, and ease of mass production. However, when compared to traditional gyroscopes, the MEMS gyroscopes are still not sufficiently accurate for many applications, such as space vehicles.
To enlarge the application field of MEMS gyroscopes, it is necessary to improve their accuracy further while not restricting their advantages (e.g., some effective modeling and processing methods were proposed in References [4,5] to reduce temperature drift). Due to the large increase in cost and time consumption of the previous methods focusing on the device itself [6,7], the gyro array was presented with the development of multisensor information fusion technology. The gyro array technique was first proposed by Bayard and Ploen [8], who integrated four gyroscope outputs to obtain an optimal rate signal by data fusion. The gyro array technique made rapid progress after that [9][10][11].
The key to the gyro array technique is the estimate method. A Kalman filter and its variants and extensions are always used, which have been proven effective to some extent. The fusion method used in the works of Bayard and Pleon was the Kalman filter. Then, an extended Kalman filter to fuse multiple gyroscopes for accuracy improvement was presented by Stubberud [12]. After that, Chang,

The Dual Uncertainties Model of a Gyro Array
The most common model for a target maneuver is the white-noise acceleration model [24]. It assumes that the target acceleration is an independent process, which is strictly white noise. It differs from the non-maneuver model only in the noise level: The white noise process w used to model the effect of the control input u has a much higher intensity than the one used in a non-maneuver model. Considering that the MEMS gyroscopes always have tiny sampling internally, this kind of model can be used to model attitude maneuver of the gyro or the gyro-installed carrier. The magnitude of the process noise is used to adjust the dynamic characteristics in the model. Due to the fact that both stochastic and set-membership uncertainties always exist simultaneously in the outputs of the gyroscopes, the angular acceleration is assumed to be an independent process with dual uncertainties. Thus, we call it a dual-noise acceleration model. Assume that ω(t) is the angular velocity of the carrier in a certain coordinate direction: The kinematic equation of the dual-noise acceleration model is then given as .
where w(t) is a zero-mean Gaussian noise (stochastic uncertainty), and d(t) is the bounded noise (set-membership uncertainty). In fact, the dual uncertainties can be interpreted as a set bounded uncertainty incorporated into the prior stochastic uncertainty [25]. For the dual-noise acceleration model given in Equation (1), the angular velocity acceleration . ω is the sum of the zero-mean Gaussian noise w ∼ N (0, σ 2 ) with variance σ 2 and the bounded uncertainty d ∈ G ⊂ R. Then . ω can be regarded as a random variable with mean d. The special thing is that it cannot be seen as a single distribution, but a set (i.e., {N (d, σ 2 )|d ∈ G }), as illustrated in Figure 1. of the gyroscopes, the angular acceleration is assumed to be an independent process with dual uncertainties. Thus, we call it a dual-noise acceleration model.
Assume that ( ) t ω is the angular velocity of the carrier in a certain coordinate direction: The kinematic equation of the dual-noise acceleration model is then given as where ( ) w t is a zero-mean Gaussian noise (stochastic uncertainty), and ( ) d t is the bounded noise (set-membership uncertainty). In fact, the dual uncertainties can be interpreted as a set bounded uncertainty incorporated into the prior stochastic uncertainty [25]. For the dual-noise acceleration model given in Equation (1) Then ω  can be regarded as a random variable with mean d . The special thing is that it cannot be seen as a single distribution, but a set (i.e.,  ), as illustrated in Figure 1. In the multidimensional case, it can be expressed as { ( , ) , } n n n n × . Then the state estimate is characterized by a set of probability densities through the combination of stochastic and set-membership uncertainties. Our proposed method poses the gyroscope filtering problem as a hidden Markov model (HMM).
Assume that the state vector is ( ) θ is the attitude angle. Considering the dual uncertainties, the discrete-time state-space model for the gyroscope array in the carrier can then be described as with the mode transition governed by a Markov chain with probabilities is the stochastic part of the process noise; In the multidimensional case, it can be expressed as {N (d, C)|d ∈ G n ⊂ R n , C ∈ R n×n }. Then the state estimate is characterized by a set of probability densities through the combination of stochastic and set-membership uncertainties.
Our proposed method poses the gyroscope filtering problem as a hidden Markov model (HMM). Assume that the state vector is x(t) = [θ(t), ω(t)] T , where θ(t) is the attitude angle. Considering the dual uncertainties, the discrete-time state-space model for the gyroscope array in the carrier can then be described as with the mode transition governed by a Markov chain with probabilities In this research, six separate gyroscopes were taken to construct a gyroscope array. According to the model shown in Equation (1), the magnitude of the process noise directly affects the ability of the model to reproduce the dynamic characteristics of the input signal. In view of this situation, this paper designed a combined filter using multiple models with different process noise magnitudes and making the noise magnitude cover a certain range by interaction: Thereby the optimality of the estimation and the stability of the filter were guaranteed. Therefore, the model set with dual uncertainties designed in this paper is given by where CorrM is a 6 × 6 correlation matrix of gyro array, r 2 is the Gaussian measurement noise variance of component gyroscopes, e 2 is a coefficient in the matrix defining the shape of the measurement noise set, and ∆T is the sampling interval. Obviously the differences between the models in Equation (7) are the values of Q i and D i , which represent the level of the stochastic part and set-membership part of the process noise.

Ellipsoidal Set and Relevant Properties
When considering the set-membership uncertainty, the set operation is involved. In order to simplify the calculations, we confined the bounding feasible sets to the ellipsoids. Several definitions and lemmas about ellipsoid sets are given below, which were used in the derivation of the multi-model combined filter with dual uncertainties.

Definition 1.
A bounded ellipsoid E of R n is defined by [23] E (a, M) = {x ∈ R n : (x − a) T where a ∈ R n is the center of E , and M ∈ R n×n is a positive-definite matrix that specifies its size and orientation. Then d k and e k in Equations (5) and (6) are actually bounded by the ellipsoids E (0, D

Definition 2.
Given an ellipsoid E (a, M)(a 1 ∈ R n ) and a full rank matrix A ∈ R n×n , then {Ax : x ∈ E (a, M)} can be expressed as AE (a, M). Lemma 1. Given an ellipsoid E (a, M)(a 1 ∈ R n ) and a matrix A ∈ R n×n , then AE (a, M) = E (Aa, AMA T ).
Proof. According to Definition 1, Assume y = Ax and then x = A −1 y, and thus After some manipulations, Equation (10) becomes This means y ∈ E A a , AMA T . The proof is complete.

Definition 3.
The Minkowski sum Ψ s of the two ellipsoids E (a 1 , M 1 ) and E (a 2 , M 2 ) is defined as denoted as The ellipsoid satisfying E (a s , M s ) ⊇ Ψ s is called the outer bounding ellipsoid of the Minkowski sum Ψ s , as illustrated in Figure 2.
denoted as ( ) ( ) The ellipsoid satisfying ( )  Proof. See the works of Maksarov and Norton [22]. For a weighted Minkowski sum of several ellipsoids, the outer bounding ellipsoid is given by the following lemma.
Proof. According to Lemma 1, it can be easily obtained that The following lemma is given to calculate the outer bounding ellipsoid of the Minkowski sum of two ellipsoids: Lemma 2. ∀p ∈ (0, +∞), E (a s , M s ) ⊇ E (a 1 , M 1 ) ⊕ E (a 2 , M 2 ), with a s = a 1 + a 2 (a 1 , a 2 ∈ R n ) and Proof. See the works of Maksarov and Norton [22].
For a weighted Minkowski sum of several ellipsoids, the outer bounding ellipsoid is given by the following lemma.

Lemma 3.
Given µ k ∈ R and the ellipsoids E (a k , M k ) (k = 1, 2, · · · r, r ∈ R + ), ∀α ∈ R r , with α k > 0 and Proof. According to Lemma 1, it can be easily obtained that On this basis and combined with relevant conclusions in Reference [23], a parametrized family of ellipsoids given by Equations (14) and (15) is obtained.

Multi-Model Combined Filtering with Dual Uncertainties (MMCF)
Following the interacting multiple model (IMM) filter [26] merging strategy, the proposed multi-model state estimate combining dual uncertainties can be updated through 4 steps: First, the mixing probabilities and the reinitialized estimates for the corresponding mode-matched filter are calculated. Then the mixing estimates are used to input into the mode-matched filter new measurements and to update the estimates. Next, the updated mode probability is calculated. Finally, the overall estimates are found using moment matching. The main difference is that the set operation is involved in the proposed algorithm, and the output contains the covariance matrix P k and an ellipsoidal set E (x k , X k ), which represents a set of means. The architecture of the MMCF algorithm is illustrated in Figure 3.
Thus, the weighted Minkowski sum can be transformed into On this basis and combined with relevant conclusions in Reference [23], a parametrized family of ellipsoids given by Equations (14) and (15) is obtained.

Multi-Model Combined Filtering with Dual Uncertainties (MMCF)
Following the interacting multiple model (IMM) filter [26] merging strategy, the proposed multimodel state estimate combining dual uncertainties can be updated through 4 steps: First, the mixing probabilities and the reinitialized estimates for the corresponding mode-matched filter are calculated. Then the mixing estimates are used to input into the mode-matched filter new measurements and to update the estimates. Next, the updated mode probability is calculated. Finally, the overall estimates are found using moment matching. The main difference is that the set operation is involved in the proposed algorithm, and the output contains the covariance matrix k P and an ellipsoidal set ( , ) k k x X  , which represents a set of means. The architecture of the MMCF algorithm is illustrated in Figure 3. Figure 3. Structure of multi-model combined filtering with dual uncertainties (MMCF) algorithm (with three models).
At time k, the MMCF algorithm can be given as follows.

Model-Conditioned Reinitialization
The predicted probability is given by where ( ) 1 j k μ − is the calculated probability for model , , , Then the mixing weight is calculated by x be the state estimate of the subfilter for model  At time k, the MMCF algorithm can be given as follows.

Model-Conditioned Reinitialization
The predicted probability is given by where µ (j) k−1 is the calculated probability for model m j k−1 at time k − 1, and Z k−1 = {z 1 , z 2 , · · · , z k−1 }. Then the mixing weight is calculated by k−1 be the error characteristics of the state estimate for the set-membership uncertainty and the stochastic uncertainty, respectively.
Suppose that only the zero-mean Gaussian uncertainty exists (i.e., X (j) k−1 = 0): Then the mixing estimate and covariance matrix are When considering the set-membership uncertainty, the mixing set of means therefore yields the Minkowski sum Then the object in this step is to find an optimal ellipsoid E ( x According to Lemma 3, there is a family of ellipsoids that meet this condition with a center x (i) k−1 and the matrix defining the shape where the parameter α k−1 ) by using the volume minimization criterion or trace minimization criterion. In this paper, the trace of X (i) k−1 was minimized due to the explicit expression of the parameter and avoidance of the nonlinear equation solution. Obviously, this criterion is favorable to computation efficiency. Then the parameter optimization problem is given by Based on the Lagrangian multiplier method, the optimal value of α

Model-Conditioned Filtering
Then each subfilter for model m i k updates state estimates with the new measurement z k . From Equations (2) and (5), the prediction set of means can be given by The ellipsoid E (x (27) and Equations (27) and (28) are obtained on the basis of Lemma 2. The parameter p that optimizes the The associated covariance matrix is stated by In the case of vanishing bounded perturbations, the state estimate is updated by means of the Kalman filter, as below: where the gain K (i) k is given by In the presence of set-membership uncertainty, the measurement z k combined with the ellipsoidal error can be regarded as an ellipsoid set, Then the predicted set E (x with the parameter q chosen to minimize the trace of X (i) k : This is

Mode Probability Update
The mode likelihood is obtained as where k . On this basis, we can calculate the model probability by (39)

Estimation Fusion
Using the matching model, the overall set of means can be expressed as Obviously the problem here is similar to that in Section 4.1, and thus the minimal-trace ellipsoid E (x k , X k ) containing the sum X k is obtained based on Equations (20)-(25): It should be noted that this solution is optimal only in the family of ellipsoids, as in Equation (23). The overall covariance is stated as Remark 1. Finally, the proposed algorithm is implemented as below: Step 1. Initializex 0 , P 0 , X 0 , and set k ← 1 .
Step 2. Model-conditioned reinitialization: Compute the mixing weight µ j|i k−1 using Equations (18) and (19), and then for each sum filter, compute x Step 4. Mode probability update: Calculate the model probability using Equation (39).
Step 5. Estimation fusion: Using the matching model, calculate the overall estimate resultsx k , X k , and P k using Equations (41)-(43), respectively.
Step 6. Set k ← k + 1 and return to Step 2.

Experimental Setup
In this section, experiments were implemented to demonstrate the performance of the proposed filter. Six MEMS gyroscope chips (ADXRS300 [27]) with external circuits were soldered onto the same printed circuit board (PCB) to construct an array. The prototype of the gyro array is shown in Figure 4. The output signals of the gyroscopes were collected through the data acquisition module constructed by a PCI (peripheral component interconnection) extension for instrumentation (PXI) system. The bandwidth of an individual gyroscope was 40 Hz, so the sampling rate was set to be 200 Hz to satisfy the Nyquist theorem.

Experimental Setup
In this section, experiments were implemented to demonstrate the performance of the proposed filter. Six MEMS gyroscope chips (ADXRS300 [27]) with external circuits were soldered onto the same printed circuit board (PCB) to construct an array. The prototype of the gyro array is shown in Figure  4. The output signals of the gyroscopes were collected through the data acquisition module constructed by a PCI (peripheral component interconnection) extension for instrumentation (PXI) system. The bandwidth of an individual gyroscope was 40 Hz, so the sampling rate was set to be 200 Hz to satisfy the Nyquist theorem.  In the experiment, the gyro array was put on a horizontal turntable. The experiment lasted for 180 s, during which the turntable was set to rotate as follows: During the other times, except for the above time periods, the turntable kept still. It should be noted that it took a period of time to adjust the motion state each time for this turntable. Thus, both the start and end phases of each motion state had transition periods, and only the middle part was in accordance with the set conditions above. The true angular rate of the turntable is shown in Figure 5. In the experiment, the gyro array was put on a horizontal turntable. The experiment lasted for 180 s, during which the turntable was set to rotate as follows: During the other times, except for the above time periods, the turntable kept still. It should be noted that it took a period of time to adjust the motion state each time for this turntable. Thus, both the start and end phases of each motion state had transition periods, and only the middle part was in accordance with the set conditions above. The true angular rate of the turntable is shown in Figure 5. Then the rate signals of the gyroscopes in the array were gathered and processed.
Through multiple tests, the correlation factor between the component gyroscopes was obtained and is shown in Table 1: It played an important role in accuracy improvement of the gyro array [28].
This meant the elements in CorrM were given. The correlation factor was analyzed with static rate signals for 1 h, and the computing method of the correlation factor followed [14].   Then the rate signals of the gyroscopes in the array were gathered and processed. Through multiple tests, the correlation factor between the component gyroscopes was obtained and is shown in Table 1: It played an important role in accuracy improvement of the gyro array [28]. This meant the elements in CorrM were given. The correlation factor was analyzed with static rate signals for 1 h, and the computing method of the correlation factor followed [14]. Four models were used in the experiment, and other parameters in the model were given as ∆T = 0.005, r = 0.6, and e = 0.6. Q i and D i for different models are shown in Table 2. Table 2. Q i and D i for different models.

Results
For comparison purposes, the stand IMM estimator using a Kalman filter and the combined filter in Section 4.2 using a single model were implemented. First, the output errors of one of the individual gyroscopes were analyzed. By using a sliding window with a width of 50, the means of the error over different time periods were calculated, as shown in Figure 6. It can be seen that the error mean value of the gyroscope changed significantly with time due to the dynamic change of input. This indicated that the assumption of two random variables as noise representation in this paper was reasonable, because a set of means was the main feature distinguishing this assumption from the single stochastic uncertainty assumption. In addition, it can be seen that the means of different time periods were all within the hypothetical set constructed by the parameters given in Section 5.1. gyroscopes were analyzed. By using a sliding window with a width of 50, the means of the error over different time periods were calculated, as shown in Figure 6. It can be seen that the error mean value of the gyroscope changed significantly with time due to the dynamic change of input. This indicated that the assumption of two random variables as noise representation in this paper was reasonable, because a set of means was the main feature distinguishing this assumption from the single stochastic uncertainty assumption. In addition, it can be seen that the means of different time periods were all within the hypothetical set constructed by the parameters given in Section 5.1. Based on this assumption, the proposed multi-model combined filter could get an ellipsoidal feasible set of means at each moment (i.e., ( , ) k k  x X ), while the standard IMM filter with respect to Gaussian noise could only obtain a point estimate. In addition, the center of the ellipsoid ˆk x could be considered to be the point estimate when needed. Figure 7 shows the state ellipsoidal estimates from 25 s to 30 s in the phase plane using the proposed algorithm. It should be noted that, for clarity, only part of the ellipsoids were selected to be shown in the figure, and the time interval between two adjacent ellipsoids was 0.05 s. The center of each ellipsoid is highlighted with "*", and the true state at the corresponding moment is highlighted with "+". It can be seen from Figure 7 that all of the values of the true state at every moment lay within the corresponding ellipsoid. On this basis, the upper and lower boundaries of the estimation were acquired by using the shape matrices of ellipsoids, as shown in Figure 8. It can be guaranteed that the angular rate estimates are within a hard boundary and that the bounded sets contain the true value, which is meaningful for the robustness of the system in control and guidance applications. Based on this assumption, the proposed multi-model combined filter could get an ellipsoidal feasible set of means at each moment (i.e., E (x k , X k )), while the standard IMM filter with respect to Gaussian noise could only obtain a point estimate. In addition, the center of the ellipsoidx k could be considered to be the point estimate when needed. Figure 7 shows the state ellipsoidal estimates from 25 s to 30 s in the phase plane using the proposed algorithm. It should be noted that, for clarity, only part of the ellipsoids were selected to be shown in the figure, and the time interval between two adjacent ellipsoids was 0.05 s. The center of each ellipsoid is highlighted with "*", and the true state at the corresponding moment is highlighted with "+". It can be seen from Figure 7 that all of the values of the true state at every moment lay within the corresponding ellipsoid. On this basis, the upper and lower boundaries of the estimation were acquired by using the shape matrices of ellipsoids, as shown in Figure 8. It can be guaranteed that the angular rate estimates are within a hard boundary and that the bounded sets contain the true value, which is meaningful for the robustness of the system in control and guidance applications. from 25 s to 30 s in the phase plane using the proposed algorithm. It should be noted that, for clarity, only part of the ellipsoids were selected to be shown in the figure, and the time interval between two adjacent ellipsoids was 0.05 s. The center of each ellipsoid is highlighted with "*", and the true state at the corresponding moment is highlighted with "+". It can be seen from Figure 7 that all of the values of the true state at every moment lay within the corresponding ellipsoid. On this basis, the upper and lower boundaries of the estimation were acquired by using the shape matrices of ellipsoids, as shown in Figure 8. It can be guaranteed that the angular rate estimates are within a hard boundary and that the bounded sets contain the true value, which is meaningful for the robustness of the system in control and guidance applications.  Below, the center of the ellipsoid was used as the point estimate, and the accuracy of the gyro array was analyzed. The original rate signals of the individual gyroscopes and processing results by different methods are given in Figure 9. The output errors and estimated errors are shown in Figure  10. In order to demonstrate the effect of the above filters clearly, the root mean square error (RMSE) of the estimated results was calculated and is presented in Table 3. For further tests, the errors were analyzed through the Allan deviation method, and the results are shown in Figure 11. Besides, the angular random walk (ARW), rate random walk (RRW), and bias instability (BI) were calculated from an Allan deviation plot, and are also presented in Table 3. In the Allan deviation test, the measurements were recorded for 1 h at a sampling rate of 200 Hz. Before that, the gyroscopes were powered on for 0.5 h. The experiments were performed in a temperature control turntable, and the temperature was set to be 25 °C. Other experimental configurations, such as hardware setup, were the same as described in Section 5.1. The turntable kept still during most times, except for 1 min. The turntable was set to swing with 10° angle amplitude and 0.5 Hz frequency in this minute to test the dynamic performance of the proposed method. Then the difference between the output of the gyroscopes (array) and the actual rotational speed of the turntable was analyzed as the output error.  Below, the center of the ellipsoid was used as the point estimate, and the accuracy of the gyro array was analyzed. The original rate signals of the individual gyroscopes and processing results by different methods are given in Figure 9. The output errors and estimated errors are shown in Figure 10. In order to demonstrate the effect of the above filters clearly, the root mean square error (RMSE) of the estimated results was calculated and is presented in Table 3. For further tests, the errors were analyzed through the Allan deviation method, and the results are shown in Figure 11. Besides, the angular random walk (ARW), rate random walk (RRW), and bias instability (BI) were calculated from an Allan deviation plot, and are also presented in Table 3. In the Allan deviation test, the measurements were recorded for 1 h at a sampling rate of 200 Hz. Before that, the gyroscopes were powered on for 0.5 h. The experiments were performed in a temperature control turntable, and the temperature was set to be 25 • C. Other experimental configurations, such as hardware setup, were the same as described in Section 5.1. The turntable kept still during most times, except for 1 min. The turntable was set to swing with 10 • angle amplitude and 0.5 Hz frequency in this minute to test the dynamic performance of the proposed method. Then the difference between the output of the gyroscopes (array) and the actual rotational speed of the turntable was analyzed as the output error.   Table 3 shows that the RMSE was reduced from 0.4916 s° to 0.1628 s° and 0.1478 s° by the IMM filter and MMCF, respectively. It reveals that the performance of the MMCF was higher than that of the IMM filter. From the Allan deviation plot and the results in Table 3, the ARW noise, RRW noise, and bias instability were all reduced by using the multi-model methods and combined filter with Models 3 and 4 to fuse measurements from the gyro array. In addition, it is obvious that the reduction factor by the MMCF was greater than that by other methods. As for the filter using Models 1 and 2, the above three kinds of noise were difficult to calculate due to large sinusoidal errors. The reason for the failure of Models 1-3 in Figure 11 was that these models could not accurately reproduce the dynamic characteristics of the input signal, and the filters using these models led to amplitude attenuation when the turntable swung. This means large sinusoidal errors were added to the output noise, similarly to the errors in Figure 10.
To further analyze the performance of the proposed method, the estimated errors over different time periods are shown in Figure 12. The detailed results are illustrated in Table 4.  In Table 3, IF is the improve factor, defined as where RMSE s refers to the RMSE of the rate signal for the single gyroscopes in the array, and RMSE a is that of the gyro array. It can be seen from these results that both the multi-model methods and the combined filter using some single models (Models 3 and 4) were effective in combining numerous gyroscopes, and the multi-model methods performed better. When using Model 1 or Model 2, the errors could be well suppressed by the filter when the input signal was 0 or constant. However, when the turntable swung, the filter led to amplitude attenuation, which demonstrated that this model could not accurately reproduce the dynamic characteristics of the input signal. When using Model 3 or Model 4, the result was the opposite: The errors were greatly reduced with sinusoidal rate input, but the filter had poor performance with constant rate input. Through switching the models according to different situations, the multi-model methods could reflect the dynamic characteristic of the input signal more accurately. Table 3 shows that the RMSE was reduced from 0.4916 • /s to 0.1628 • /s and 0.1478 • /s by the IMM filter and MMCF, respectively. It reveals that the performance of the MMCF was higher than that of the IMM filter. From the Allan deviation plot and the results in Table 3, the ARW noise, RRW noise, and bias instability were all reduced by using the multi-model methods and combined filter with Models 3 and 4 to fuse measurements from the gyro array. In addition, it is obvious that the reduction factor by the MMCF was greater than that by other methods. As for the filter using Models 1 and 2, the above three kinds of noise were difficult to calculate due to large sinusoidal errors. The reason for the failure of Models 1-3 in Figure 11 was that these models could not accurately reproduce the dynamic characteristics of the input signal, and the filters using these models led to amplitude attenuation when the turntable swung. This means large sinusoidal errors were added to the output noise, similarly to the errors in Figure 10.
To further analyze the performance of the proposed method, the estimated errors over different time periods are shown in Figure 12. The detailed results are illustrated in Table 4. Models 3 and 4 to fuse measurements from the gyro array. In addition, it is obvious that the reduction factor by the MMCF was greater than that by other methods. As for the filter using Models 1 and 2, the above three kinds of noise were difficult to calculate due to large sinusoidal errors. The reason for the failure of Models 1-3 in Figure 11 was that these models could not accurately reproduce the dynamic characteristics of the input signal, and the filters using these models led to amplitude attenuation when the turntable swung. This means large sinusoidal errors were added to the output noise, similarly to the errors in Figure 10.
To further analyze the performance of the proposed method, the estimated errors over different time periods are shown in Figure 12. The detailed results are illustrated in Table 4.   Table 4 also demonstrate the performance improvements.
It can be seen from Table 4 that the RMSEs of the estimation were obviously reduced compared to the original rate signals of the individual gyroscopes. Furthermore, the performance of the proposed method had to do with the dynamic characteristics of the input signal. The reduction factor with respect to 0 input and constant input is about 4-5, whereas the corresponding factor for sinusoidal input was about 2-3, which was lower than that of the 0 input and constant input. Additionally, it indicated an increase in the RMSE with increasing frequency and amplitude with respect to the swing rate signal.    Figure 12 shows that the proposed method could effectively reduce the errors of the gyroscope in any kind of motion. The statistics results in Table 4 also demonstrate the performance improvements. It can be seen from Table 4 that the RMSEs of the estimation were obviously reduced compared to the original rate signals of the individual gyroscopes. Furthermore, the performance of the proposed method had to do with the dynamic characteristics of the input signal. The reduction factor with respect to 0 input and constant input is about 4-5, whereas the corresponding factor for sinusoidal input was about 2-3, which was lower than that of the 0 input and constant input. Additionally, it indicated an increase in the RMSE with increasing frequency and amplitude with respect to the swing rate signal.

Conclusions
In this paper, a multi-model combined filter allowing for simultaneous treatment of stochastic and set-membership uncertainties was proposed to combine multiple MEMS gyroscopes to improve overall accuracy. The dual uncertainties model of a gyroscope array was established. To reproduce the dynamic characteristics of the rate signal, multiple models with different process noises were involved in the filter. Combined with the IMM update process and the calculation method for a Minkowski sum of multiple ellipsoidal sets, the multi-model state estimates integrating dual uncertainties was updated.
The experimental results showed that the RMSE of the estimated rate signal by the proposed method was reduced from 0.4916 • /s to 0.1478 • /s, which performed better than the IMM filter and combined filters with single models. This proved that the proposed method is efficient in improving the system overall performance and that it produces better adaptability to different kinds of uncertainties and different dynamic characteristics. Moreover, the proposed method can provide the bounds of the rate signal estimates. This property is meaningful for the robustness of the system and benefits attitude control and guidance.