Reduced-Drift Virtual Gyro from an Array of Low-Cost Gyros

A Kalman filter approach for combining the outputs of an array of high-drift gyros to obtain a virtual lower-drift gyro has been known in the literature for more than a decade. The success of this approach depends on the correlations of the random drift components of the individual gyros. However, no method of estimating these correlations has appeared in the literature. This paper presents an algorithm for obtaining the statistical model for an array of gyros, including the cross-correlations of the individual random drift components. In order to obtain this model, a new statistic, called the “Allan covariance” between two gyros, is introduced. The gyro array model can be used to obtain the Kalman filter-based (KFB) virtual gyro. Instead, we consider a virtual gyro obtained by taking a linear combination of individual gyro outputs. The gyro array model is used to calculate the optimal coefficients, as well as to derive a formula for the drift of the resulting virtual gyro. The drift formula for the optimal linear combination (OLC) virtual gyro is identical to that previously derived for the KFB virtual gyro. Thus, a Kalman filter is not necessary to obtain a minimum drift virtual gyro. The theoretical results of this paper are demonstrated using simulated as well as experimental data. In experimental results with a 28-gyro array, the OLC virtual gyro has a drift spectral density 40 times smaller than that obtained by taking the average of the gyro signals.


Introduction
The concept of combining the outputs of an array of high-drift sensors to produce a virtual low-drift sensor was introduced by Bayard and Ploen [1,2], who used a Kalman filter to produce the virtual sensor. They showed, through analysis and simulation examples, that the reduced drift of the virtual sensor depends on the drift components of the individual sensors being "favorably correlated" (i.e., having some negative correlation). In order to implement the Bayard approach it is necessary to estimate these correlations from sensor data, but no such estimation algorithm has appeared in the literature. Several researchers have attempted to implement the Bayard approach (see [3][4][5][6]). However, none of these efforts were able to implement the full Bayard algorithm because they were not able to estimate the correlations between the drift components of the individual sensors. In this paper we present an algorithm for estimating these correlations.
In addition to providing the Kalman filter algorithm for combining sensors, Bayard and Ploen also derived a mathematical formula for the drift of their optimal virtual sensor. This formula provides a lower bound on the drift of any method of combining the sensor outputs. In this paper we consider obtaining a virtual sensor by taking a linear combination of the individual sensor outputs. We derive a formula for computing the optimal coefficients and also derive a formula for the drift of the resulting virtual sensor. Our formula for the drift of the optimal linear combination virtual sensor is identical to Bayard and Ploen's formula for the drift of the Kalman filter-based virtual sensor. Because the optimal drift can be obtained using only a linear combination of sensor outputs, we conclude that a Kalman filter is unnecessary for this problem.
The performance of a rate gyro or accelerometer is influenced by a number of factors including scale-factor errors, quantization effects, temperature effects, random drift, and additive noise. A comprehensive account of all of these factors for gyros is found in an IEEE Standard [7]. Among the many possible noise sources mentioned, the additive noise and random drift are often the dominant sources. In our experience, this is the case for tactical-grade MEMS sensors. Other researchers have made the same observation [8].
When an array of multiple sensors is considered, there may be correlations in the statistical properties of one sensor to another that can be exploited to reduce the drift. In this paper, we assume that the additive noise terms for different sensors are uncorrelated but that there may be correlations between the drift components of different sensors. This paper provides an algorithm to estimate these correlations from a finite amount of calibration data. We consider an array of gyros, but the modeling algorithm applies also to an array of accelerometers [9,10]. The organization of this paper, as well as its relationship to previous work, is as follows.
In Section 2.1 the standard single-gyro model for additive noise (ARW) and random drift (RRW) is introduced, along with the concept of Allan variance, which is computed in discrete-time using sampled gyro signals. The approach described here is based on the computed Allan variance of a single gyro, but is not based on the usual method of fitting lines to the Allan variance plot [7,11]. It was shown in [12] that this line-fitting method is a statistically inconsistent estimator. In [13] a modification of Allan variance using sliding windows was presented in order to produce a better plot. However, this modification is not necessary for the statistical modeling procedure used here because the procedure automatically deemphasizes inaccurate Allan variance points using theoretical formulas for the variance and covariance of computed Allan variance points. Other approaches to modeling individual gyros are presented in [14,15]. However, these approaches have not been extended to modeling the correlations in an array of gyros. In [16], the formulas of [10] were used to model triaxial rate gyros. However, the three gyro signals were assumed to be uncorrelated. Section 2.4 presents a model for the signals in a two-gyro array, and introduces a new statistic that we call the Allan covariance between two gyros. The model for the correlation between two gyros can be extended to an array of many gyros by modeling all pairwise combinations of gyros. The statistical properties of a gyro array are described by the spectral density matrices, R and Q of the additive noise and random drift components of the gyro signals. Section 2.7 extends to 2-gyro model to an array of gyros and derives an algorithm for estimating the off-diagonal elements of the Q matrix. Also derived is a formula for the coefficient vector c * , which gives the optimal linear combination of gyro outputs to produce a scalar virtual gyro signal. Section 2.9 presents simulation results for a 6-gyro array. Signals are generated using given Q and R matrices. EstimatesQ andR are obtained from a finite amount of data using the algorithm from Section IV. The matrixQ matrix is used to calculate the coefficients for the optimal linear combination of gyro signals. The theoretical formula for the drift of the virtual gyro is confirmed by simulations. Section 2.11 presents experimental results with an array of 28 MEMS gyros. The algorithm given in Section 2.7 is used to estimate the Q and R matrices for this gyro array and to compute the optimal coefficient vector c * . The RRW spectral density of the resulting virtual gyro is calculated and found to be over 40-times smaller than that for a virtual gyro obtained by averaging the gyro outputs.

Single-Gyro Model for Noise and Drift
The output signal, y(t), of a gyro in the absence of motion will be referred to as a calibration signal, and can be described by the following equations: where b(t) is the bias (random drift), n(t) is observation noise, which is assumed to be white noise with spectral density R, and w(t) is another white-noise process described below in Equation (2). In the discrete-time treatment given below it is assumed that the gyro signal is lowpass filtered to a bandwidth B Hz. For rate gyros, the bias term b(t) is referred to as rate random walk or RRW, and the additive noise term n(t) is referred to as angle random walk or ARW. The bias term b(t) is a random process consisting of the integral of a white noise process w(t) with spectral density Q. The integral creates a nonstationary random process whose variance increases with time, which provides a good model for the drift of a sensor. The bias signal is modeled as followṡ The bias and noise terms are characterized by the numbers Q and R, respectively. Allan variance, denoted a(τ), is computed for different values of a smoothing parameter, τ. In order to compute the Allan variance associated with a record of calibration data, the signal is partitioned into N τ-second intervals.
Letȳ k , k = 1, · · · , N be the average values of the N nonoverlapping τ-second intervals of calibration data. Then the Allan variance statistic for this data set is defined as The sensor signal is sampled with a sampling interval of T seconds to obtain the discrete-time signal Consider the statistical characterization of the sampled noise and drift components of the sensor signal. The noise term n(t) is modeled as bandlimited white noise with spectral density R and bandwidth B Hz. If the noise term is sampled with sampling interval T = 1/(2B) second, the noise samples have variance R/T and noise samples at different sampling instants are uncorrelated. That is, where δ[k] = 1 if k = 0 and zero otherwise. The drift term, b(t), is modeled as integrated white noise, as shown in Equation (2). The sampled bias term may be written as: Let the random variables z k , k = 1, · · · , M, be the averages of M nonoverlapping length-m segments of data, where M = N/m: In this analysis, the Allan variance smoothing parameter, τ, is τ = mT.
The Allan variance statistic is computed from a finite record of sampled data as follows:

Mean and Covariance
The results in this subsection were derived in [10]. The mean values for the ARW and RRW components of a[m] are: Consider now the Allan variances computed from a data record of N samples at two values of m, say m 1 and m 2 , where m 2 = p · m 1 for some integer p ≥ 1. We require that M 1 and M 2 defined as In order to use Equations (11) and (12) for a set of Allan variance points, the values of m used to compute these points must all be integer multiples of each other. This is accomplished by using values of m that are powers of 2. Consider a set of points m = 2 4 · · · 2 q for some integer q. Let a[m] be the corresponding vector of Allan variance statistics. The covariance matrices of a[m] for the ARW and RRW terms are obtained, element by element, from Equations (11) and (12), respectively. We denote these covariance matrices as covar(a[m]) = C R (m, R) for ARW (13) covar(a[m]) = C Q (m, Q) for RRW.
For example, the ij element of the matrix C R (m, R) is covar(a[m i ], a[m j ]), where m i and m j are the ith and jth elements, respectively of m and Equation (11) is used to evaluate the covariance.

Algorithm for Estimating Q and R for a Single Gyro
The algorithm presented in this section is based on standard results in statistical estimation theory, as expressed by the following theorem describing the best linear unbiased estimator [17]. This estimator has the lowest variance of any linear unbiased estimator.
where x is an N × 1 vector of observations, H is a known N × p matrix, θ is a p × 1 vector of parameters to be estimated, and w is zero-mean random vector with covariance matrix C. Then the best linear unbiased estimator The expected value of the estimator is the true parameter vector, θ, and the covariance matrix of the estimator is This theorem was used in [10] to obtain an algorithm for estimating Q and R for a single gyro and is used in Section IV below to derive an algorithm for modeling an array of gyros. The single-gyro algorithm is summarized in Algorithm 1.

Algorithm 1:
An algorithm for estimating individual gyro parameters Q and R. c 2012 IEEE, reprinted with permission from [10]. 1. Given a data record of N points of gyro calibration data let J = J 1 − 3, where J 1 is the integer part of log 2 (N). Calculate the Allan variance statistic for m = 2 4 · · · 2 J T .
2. Letτ 0 = Tm 0 be an estimate of the smoothing lag corresponding to the minimum value of the Allan variance. For example, choose m 0 to be the index from Step 1 that gives the smallest Allan variance.

From the values of m calculated in
Step 1, select m 1 = 2 · · · 2 N 1 T , such that the maximum value 5. Calculate a preliminary estimate of Q:Q 0 = 3R 0 /τ 2 0 . 6. Let (see Equations (13) and (14)): where m is obtained in Step 1. 7. The estimates of Q and R are given by

Mathematical Model for a Two-Gyro Array
In order to model an array of gyros it suffices to jointly model two gyros by estimating the correlation properties between the two gyros, as well as their individual statistical properties. Once the procedure is developed for two gyros it can be used for an array of any number of gyros by modeling all pairwise combinations of gyros.
The output signal vector, y(t), of a two-gyro array can be described by the following equation where ω(t) is a vector of true rate signals, b(t) is a vector of gyro biases, and n(t) is an observation noise vector. In the absence of motion, the true rate is zero and the gyro output signal consists of bias plus noise. The bias term b(t) is a vector of rate random walks or RRW, and the additive noise term n(t) is a vector of angle random walks or ARW.
The elements of the noise vector n(t) are modeled as statistically independent bandlimited white noise processes with spectral densities R 1 and R 2 , respectively, and common bandwidth B Hz. The spectral density matrix of n(t) is Let w(t) be a vector white-noise process with spectral density matrix Q. That is The bias vector is a random walk arising from w(t): Gyro calibration data, with ω(t) = 0, consists of bias and noise terms. These terms are characterized by the matrices Q and R, respectively. Given a finite record of gyro calibration data, the problem at hand is to estimate the nonzero elements of Q and R. In Algorithm 1 of the previous section, an algorithm is given for modeling an individual gyro. This algorithm may be used separately for gyro 1 and gyro 2 to obtain estimates of Q 11 and R 1 from gyro 1 data, and Q 22 and R 2 from gyro 2 data. All that remains is to show how to estimate Q 12 using the data from gyros 1 and 2.
Consider uniformly sampling the two gyros with sampling interval T: Similarly to Equation (7) for a single gyro, define the vectors z k to be the average of M nonoverlapping length-m segments of vector-value gyro data y k , where M = N/m: We define the Allan covariance matrix for a smoothing interval of m samples to be where M = N/m. Equation (23) is a simple extension of the single-gyro Allan variance formula (8). To our knowledge, this paper is the first time Equation (23) has been used to model an array of correlated gyros. The diagonal elements of A[m] are the usual Allan variances for gyros 1 and 2, respectively. The off diagonal terms, which are equal, are the Allan covariance between gyro 1 and gyro 2. We denote the off-diagonal term of A[m] as c[m], which is given by where z k (1) is the first element of z k containing signals from gyro 1, and z k (2) is the second element of z k containing signals from gyro 2. The terms in the summation defining c[m] may be written using the following vectors of signal differences, for k = 1, · · · , M − 1: . . .
We can write where o is a vector of ones. Similarly, we can write Then

Mean of the Allan Covariance
The mean of c[m] is found by calculating the expected value of x k w T k in Equation (28) The lth elements of the vectors x k and w k are given by Using these expressions, the ij element of the matrixR kk is calculated as follows: The τ 1 and τ 2 integrals are each computed over intervals of length mT. The value of the double integral is computed by counting the number of length-T overlapping intervals in the τ 1 and τ 2 integrals. When i = j, the limits of integration overlap completely, giving a value of mT for the double integral. When j = i + 1 or i − 1 there are m − 1 overlapping intervals and the value of the double integral is (m − 1)T. In general, the number of overlapping intervals is m − |j − i|. Thus, for the RRW component,R kk is equal to Using Equation (28), the mean of c[m] is The term in the summation is simply the sum of all the elements of the matrixR kk , which can be calculated using Equation (32) to be In practical applications the RRW term is significant only for large value of the smoothing parameter τ = mT, e.g., m > 1000. Thus, only the highest power of m in Equation (34) Because the mean of the ARW component of Allan variance is zero, the previous equation is also the mean of the Allan covariance of a gyro calibration signal.

Covariance of the Allan Covariance
The covariance between two values of the Allan covariance statistic, c[m 1 ] and c[m 2 ], is needed to obtain the best linear unbiased estimate of Q 12 using an algorithm similar to that used for modeling a single gyro. Let m 2 = p · m 1 for some integer p ≥ 1. Let N be the number of samples in the sampled gyro calibration signals and let M 1 and M 2 be defined as By choosing m 1 , m 2 , and N to be powers of 2, both M 1 and M 2 will be integers. The covariance expressions for the ARW and RRW components are derived in the Appendix A. For the ARW component, the covariance is Given a vector m containing n values of m (each a power of 2), let the vector a be the Allan covariance evaluated at each element of m. Then the n × n covariance matrix of a is calculated element-by-element using Equation (37 The n × n covariance matrix of a is calculated element-by-element using Equation (38). The resulting covariance matrix is called C Q (m, Q 11 , Q 22 , Q 12 ).

Statistical Model for a Gyro Array
Consider an array of g gyros. The Allan covariance matrix A[m] will be a g × g symmetric matrix. The off-diagonal terms, A ij [m] i = j are the Allan covariances between gyro i and gyro j. We assume that the ARW components of the different gyros are uncorrelated so that the off-diagonal terms of the spectral density matrix R are equal to zero. This uncorrelated assumption has been verified for the hardware gyro described in Section 2.11. Each of the diagonal elements of R may be estimated using the single-gyro algorithm given in Algorithm 1, which also gives the diagonal elements of the spectral density matrix Q. All that remains is to estimate the off-diagonal terms of Q, Q ij i = j. An algorithm to do this is derived next.
The algorithm for estimating Q ij is based on Theorem 1. In order to use this theorem, let m = 2 4 · · · 2 J be a set of J smoothing indices and let a be the vector of corresponding Allan covariances From Equation (34) the mean of a is The covariance of a is (see Equations (37) and (38)) Knowing the mean and covariance of a, Theorem 1 may be used to obtain an estimator for Q ij as follows. Using the notation of Theorem 1, let x = a, H = mT/3, θ = Q ij , with C given by Equation (41).
The only difficulty is that the matrix C Q in Equation (41) depends on Q ij , the parameter we want to estimate. The simplest way to obtain a realizable algorithm is to set Q ij to zero in the formula for C Q . The effect of this approximation will be examined in the simulation study in Section 2.9. The complete algorithm for estimating the off-diagonal entries of the matrix Q is given in Algorithm 2.

Algorithm 2:
An algorithm for estimating the off-diagonal terms, Q ij , of the spectral density matrix Q for the drift components of an array of gyros.
1. Given a data record of N points of vector-valued gyro calibration data from an array of g gyros, let J = J 1 − 3, where J 1 is the integer part of log 2 (N). Calculate the set of g × g Allan covariance matrices {A[m]} using Equations (21) through (23) for all m ∈ m, where m = 2 4 · · · 2 J T .

Calculate the covariance matrix (see Equations
6. The estimate of Q ij is calculated as follows:

Optimal Linear Combination of Gyro Signals
Suppose we have an array of g gyros with a known RRW matrix Q. The outputs of the gyros are collected into a vector y(t). A virtual gyro signal v(t) can be obtained by taking a fixed linear combination of the gyro signals: where the coefficient vector Each gyro is measuring the same velocity signal ω(t). The constraint in Equation (43) implies that the weighted signal components in each of the gyro signals will add up to ω(t).
The signal v(t) can be modeled as the output of a single gyro with some effective RRW spectral density matrix Q v . We define the optimal coefficient vector c * to be the one that results in the virtual gyro having the smallest value of Q v .
From Equations (20) and (42), the bias signal for the virtual gyro satisfieṡ Comparing Equation (44) with Equation (2), it can be seen that the spectral density of the RRW component of the virtual gyro is Thus, Q v for a virtual gyro signal formed as a fixed linear combination of gyro signals is The optimal coefficient vector c * is the one that minimizes the drift of the virtual gyro Equation (46) subject to the constraint in Equation (43). That is, Using a Lagrange multiplier to convert this to an unconstrained optimization problem and setting the first derivative to zero yields the well-known solution to Equation (47): The virtual gyro signal obtained by taking the linear combination of gyro signals using the coefficient vector c * will have an RRW spectral density Q * v given by Equation (46). Substituting Equation (48) in Equation (46) yields Because o is a vector of ones, the previous equation says that Q * v , the RRW spectral density of the optimal linear combination (OLC) virtual gyro, is equal to the reciprocal of the sum-of-elements of the matrix Q −1 . This is identical to the virtual gyro drift obtained from the Kalman filter approach. In column 10 of Bayard and Ploen's patent [2] they say that "the theoretically minimum drift attainable from combining N gyro devices having Rate Random Walk matrix Q is given by 1/(sum of the elements of Q −1 )." Thus the minimum drift achieved by the Kalman filter based (KFB) approach is identical to that achieved by the OLC approach presented in this paper.
Consider the following similarities and differences between the KFB and OLC approaches. In both cases the algorithms are derived for gyros that are modeled using only ARW and RRW components. In both cases the matrix Q must be estimated from calibration data that is obtained while the array is motionless. In order to solve the Kalman filter equations the estimated matrixQ must be positive definite. For the OLC approach, the theoretical formulas assume a positive definiteQ but it is shown in Equation (54) (and preceding text) how to handle the case whenQ is not positive definite. Finally, to calculate the KFB virtual gyro output from an array of g gyros in an operational environment, the gyro signals must be processed by an gth-order Kalman filter. In contrast the OLC virtual gyro simply forms a weighted sum of the g gyro signals using a fixed linear combination.

Simulation Example
The example in this section consists of a simulation of an array of six gyros. The matrix R is a diagonal matrix with units of deg 2 /(hr) and the following diagonal elements: We begin with some theoretical calculations based on the known Q matrix. Then, in the second subsection below, we present results using simulated gyro signals.

Theoretical Calculations
The drift quality of the individual gyros is given by the diagonal elements of Q. Gyro 3 is the worst, with Q 33 = 163 × 10 −3 . Gyro 1 is the best, with Q 11 = 11.9 × 10 −3 .
There are several ways that a virtual gyro signal can be obtained by taking a fixed linear combination of the outputs of the six gyros. One way is to take the average of the gyro signals using the coefficient vector c = 1/6 · · · 1/6 T . Another way is to take a weighted average of the gyro signals in which the weights are the inverses of the Q values for each individual gyro. This is equivalent to using only the diagonal elements of the matrix Q in Equation (48). A third way is to use the optimal linear combination c * given by Equation (48). The RRW spectral density Q v for any virtual array created by taking a fixed linear combination of gyro signals can be computed using Equation (46). The results for the three methods just mentioned are shown in Table 1. Table 1. Coefficient vectors for a 6-gyro array and corresponding RRW spectral density Q v using three methods. Method 1 is the unweighted average of the gyro signals. Method 2 uses only the diagonal elements of the theoretical Q matrix to calculate the coefficients. Method 3 is the optimal linear combination, which is calculated using the entire Q matrix.

Method 1 Method 2 Method 3
Coefficient Vector Several observations can be made about this example. First, the simple average of the gyro signals does not improve the drift of the best individual gyro, mostly because the individual gyros have very different Q values. This disparity in Q values is accounted for using only the diagonal elements of Q in the calculation of the coefficient vector, which greatly improves the value of Q v . Finally, the optimal linear combination, which is calculated using the entire Q matrix, gives the lowest value of Q v .

Simulation Results
In order to test Algorithm 2, 500 realizations of the gyro signals were created at 10 Hz sampling rate for a duration of 31.1 hours using the spectral density matrices given in Equations (50) and (51). For each of the 500 trials the matricesQ andR were computed using Algorithm 2. Two different coefficient vectors were computed. The first used the diagonal elements ofQ in Equation (48). The second, optimal coefficient vector, was computed using the fullQ matrix in Equation (48). Finally, virtual gyro signals were created by taking linear combinations of the six gyro signals using the two coefficient vectors. The virtual gyro signals were modeled using Algorithm 1 and the spectral densities (R v ,Q v ) for the virtual gyros were recorded. Figure 1 is a plot of theQ v values corresponding to the two coefficient vectors for each trial. For the virtual gyros created using only the diagonal elements of Q, the resultingQ v values had a mean of 3.9 × 10 −3 and a standard deviation of 2.9 × 10 −4 . The mean ofQ v is 3% larger than the theoretical value shown in Table 1 under Method 2. TheQ v values for the virtual gyros created using the fullQ matrix had a mean of 3.0 × 10 −3 and a standard deviation of 2.5 × 10 −4 . This mean value is 11% larger than the theoretical value shown in Table 1 under Method 3 due to errors in estimating the the off-diagonal elements of Q. We believe that those errors are due primarily to the approximation used in computing C Q . See Equation (41) and the discussion following. Coefficient vector calculated using entire Qhat matrix Coefficent vector calculated using only diagonal elements of Qhat Figure 1. Random drift (RRW) spectral densityQ v for virtual gyros corresponding to coefficient vectors calculated using only the diagonal elements ofQ and using the entireQ matrix. A virtual gyro created by averaging the gyro signals has an RRW spectral density of 11.5 × 10 −3 for all trials (see Table 1).

Hardware Results for a 28-Gyro Array
The algorithm has been implemented on a sensor array developed by NUWC. The sensor array consisted of 28 3-axis MEMS gyroscopes mounted on a flat board. The gyroscope part number is L3GD20 from ST Microelectonics. Calibration data were collected while the gyro array was motionless (lying on a table) at room temperature. The outputs of each gyro were lowpass filtered and sampled at 10 Hz. Each pair of 14 sensors are handled by a Propeller multicore microcontroller. The XYZ data from each all the sensors are packaged, time tagged and streamed to the computer through a high-speed USB port in binary format. Dedicated software on the PC was used to log the data and convert to ASCII upon saving it to a file.
The gyro signals were recorded for 58.3 hours and the x-axis signals were analyzed. In order to test the assumption that the ARW components of different gyros are uncorrelated, the Algorithm 2 was modified to estimate the entireR, including the off-diagonal elements, as well as the matrixQ. The diagonal elements ofR were sorted and compared with the eigenvalues ofR, which were also sorted. The entries in the two lists varied only in the third significant figures, indicating that the ARW components are indeed uncorrelated. The remainder of this section uses the Algorithm 2, which assumes that R is a diagonal matrix. By not explicitly estimating the off-diagonal terms in R there are fewer parameters to estimate, which results in a more accurate estimate of the off-diagonal elements of Q.
The sorted eigenvalues ofQ and the sorted diagonal elements ofQ are plotted in Figure 2. The differences in the elements of these two lists indicates correlation between the RRW components of the gyros, which can be used for drift reduction. However, Figure 2 also reveals a difficulty that must be addressed before computing the optimal linear combination vector c * fromQ. The difficulty is thatQ has several negative eigenvalues, which are due to errors in estimating the individual entries in theQ matrix. IfQ has negative eigenvalues then the quantity c TQ c, which is minimized to obtain the optimal linear combination coefficient vector, may take negative values. In this case, the result of the minimization will be the largest possible negative value of c TQ c, which is not a useful result. There are several ways to deal with this difficulty, and we have found the following to be the most effective. Let the singular value decomposition (SVD) ofQ be written as follows: where u k and v k are the left and right singular vectors, respectively, ofQ, and s k are the singular values. The inverse ofQ may be written as:Q WhenQ has negative eigenvalues, we propose to compute a partial inverse by omitting the first term or two in the SVD expansion ofQ −1 . These are the least important terms because they contain the inverse of the largest singular values. For theQ obtained from the 28-gyro array only the first term was omitted, and the coefficients for the optimal linear combination were computed as follows, using x as a temporary variable (see Equation (48)): A virtual gyro signal was created by taking a linear combination of the 28 gyro signals using three different methods. Method 1 is the unweighted average of the gyro signals. Method 2 uses only the diagonal elements ofQ in Equation (48) to calculate the coefficients. Method 3 is the optimal linear combination, which is calculated using Equation (54). The Allan variance plots for these three virtual gyros are shown in Figure 3 and the estimated R and Q parameters are shown in Table 2.    Table 2 shows that, for this set of experimental data, using only the diagonal elements ofQ gives a virtual gyro with a value ofQ v that is a factor of 8-times smaller than that obtained by an unweighted average of the gyro signals. Using the fullQ matrix (with partial inverse) decreasesQ v by an additional factor of 5.
Note that if the RRW components of the 28 gyros in the array were uncorrelated, the theoretical matrix Q would be a diagonal matrix. Due to estimation errors, the matrixQ would contain nonzero off-diagonal terms, but these off-diagonal terms would not contain any useful information. The fact that the linear combination obtained using the fullQ matrix gave substantially better results than the linear combination obtained using only the diagonal elements ofQ is evidence that the off-diagonal terms of the theoretical Q matrix are nonzero. That is, there is evidence that the RRW components for the gyros in this hardware array are correlated.

Discussion and Conclusions
We have presented an algorithm for estimating R and Q, the spectral density matrices for the ARW and RRW components, respectively, of an array of gyros. The algorithm is based on the definition of the Allan covariance between two gyros, which is a generalization of the Allan variance for a single gyro. We considered virtual gyros, defined by fixed linear combinations of the signals from a gyro array, and derived a formula for the coefficients of the optimal linear combination resulting in the virtual gyro with smallest drift.
The performance of the minimum-drift virtual gyro was demonstrated using simulated as well as experimental data. Experimental results with an array of 28 gyros showed that the proposed virtual gyro had a drift spectral density 40-times smaller than that obtained by taking the average of the gyro signals.
One area of future work is to extend the approach presented in this paper to handle gyros who calibration signals are modeled by more than just ARW and RRW components. For example, references [14,15] show how to model a single gyro which has ARW, RRW, and Gauss-Markov components. Although all of these components are needed to accurately model the gyro, the long-term drift of a gyro is determined only by the RRW component. If the modeling procedure in [14,15], or some other modeling procedure, were extended to a gyro array, only the spectral density matrix Q of the RRW components would be needed to obtain a minimum-drift OCL virtual gyro. The optimal coefficients would still be calculated using Equation (48).
In order for the minimum-drift gyro to be much better than the average of the gyros, the off-diagonal entries in the matrix Q must be nonzero; that is, there must be correlations among the RRW components of different gyros in the array. Another important area of future work is to determine what aspects of gyro fabrication influence the correlation properties between gyros. Finally, we mention that other approaches to combining multiple gyros and accelerometers are presented in [18]. It is possible that the statistical analysis presented above could be combined with the ideas in [18].
Author Contributions: Richard Vaccaro developed the mathematical analysis and derived the algorithms. Ahmed Zaki conceived of an alternative to a Kalman filter for virtual gyro synthesis and developed the hardware and software used to collect the experimental data.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Derivation of Covariance of Allan Covariance Between Two Gyros
Consider the second moment of the Allan covariance c[m] defined in Equation (28).
The expectation inside the double summation contains the product of four Gaussian random variables. Consider the following formula for such an expectation [19,20]: The correspondence between Equations (A1) and (A2) is as follows: Each of the random variables in Equation (A3) has zero mean, thus the last term of Equation (A2) is equal to zero. The first term of Equation (A2) , using Equations (A3) and (29), is Substituting this into the right-hand side of Equation (A1) and using Equation (33)  . These two terms will be called "term 2" and "term 3," respectively. They are given by and term 3 = 1 and R 1 ij and R 2 ij are defined as follows: The variance of c[m] is the sum of terms 2 and 3, and the values of these terms are computed next.
The following single-gyro formulas for the variance of the Allan variance a[m] (see Equation (8)) were derived in [10].
var(a[m]) = 1 where R ij = E(x i x j ). (A13) Term 3 in Equation (A6) for ARW is zero because additive noises on different gyros are uncorrelated. In order to compute term 3 for RRW, note that the sum of elements ofR ij is equal to the sum of elements ofR T ij . Thus, the summand of Equation (A6) may be written as (o TR ij o) 2 . Now compare the definition ofR ij in Equation (A7) with the definition of R ij in (A10). In Equation (A10) the expectation is taken of a product of signal difference vectors from a single gyro. When the data in the signal difference vectors come from overlapping intervals, the expectation is Q times the length of the overlap. In Equation (A7) the expectation is taken of the product of signal difference vectors from two different gyros. When the data in the signal difference vectors come from overlapping intervals the expectation is Q 12 times the length of the overlap. The intervals in Equation (A9) are the same lengths as those in Equation (A7). The expected values are squared in Equation (A9) as well as in Equation (A6). Thus, term 3 for RRW is given by the variance formula (A12) with the following changes: Q 2 is replaced by ARW: R 2 is replaced by R 1 R 2 /2 RRW: Q 2 is replaced by (Q 1 Q 2 + Q 2 12 )/2. , where m 2 = p · m 1 for some integer p ≥ 1. The covariance calculations are similar to the variance calculations given above except that m is replaced by m 1 for the first gyro and m 2 for the second gyro. In [10] the same intervals of length m 1 and m 2 were used with data from a single gyro. As was seen in the variance calculations, the formulas for a single gyro may be used to get the results for two gyros using the substitutions shown in Equation (A15). Because the lengths of the overlapping intervals depend on m 2 = p · m 1 and are the same for the single gyro and two gyro cases, the expressions for covar(c[m 1 ], c[p · m 1 ]) may be obtained from Equations (11) and (12)