High accuracy inertial sensors from inexpensive components

Abstract

An analytical apparatus for optimally combining measurements from N individual rate sensing devices or gyros into a single rate estimate significantly improves performance over that of any individual component device. Kalman filtering is used to combine rate sensed devices optimally in the sense of minimizing the variance of the rate error. The Riccati differential equation (RDE) associated with combining a collection of rate sensed devices is completely and exactly solved to derive to the matrix RDE. This analytic solution serves as the key for understanding all of the theoretical properties of the optimal filter, and provides a complete characterization of the final virtual rate sensed performance. In addition, the analytic RDE solution allows many practical problems to be solved that have proved essential for developing successful filter implementations. A discrete-time minimum variance filter implementation combines sensor measurements optimally.

Description

BACKGROUND OF THE INVENTION

[0002] 1. Field of the Invention

[0003] The invention relates to the field of gyroscopic instrumentation and in particular to processing the information to improve accuracy of the gyroscopic data.

[0004] 2. Description of the Prior Art

[0005] Small, lightweight, and low-power devices are needed for accurate inertial sensing in a wide range of planned future NASA missions requiring guidance, control and inertial navigation functions (e.g., Mars missions, small body landing, rendezvous and docking, planetary ascent, precision landing, atmospheric entry, aeromaneuvering, aerocapture, etc.). and reduce mass and power requirements in emerging low-cost microspacecraft missions.

[0006] Emerging microelectromechanical systems (MEMS) are known to have smaller size, lower weight and less power consumption than their discrete counterparts, and in many ways hold the potential for providing the desired capabilities. Unfortunately, present state-of-the-art MEMS devices have low-grade performance and cannot compete with established sensors in such high-accuracy application areas. However, one important aspect of MEMS devices is that they are small and can be fabricated in large numbers on a single chip. One object of the invention is to create a “virtual sensor” by combining many individual inertial sensing devices (e.g., gyros/accelerometers) into a single effective device. The resulting device would have performance properties surpassing the properties of any individual device.

[0007] Allan, Ashby and Hodge attempted in 1998 with some success to combine three inexpensive Sports stopwatches D. W. Allan et.al., “Fine-Tuning Time In The Space Age,” IEEE Spectrum, pp. 42-51, March 1998.

BRIEF SUMMARY OF THE INVENTION

[0008] The invention is a systematic method to combine multiple gyro devices so that the collective performance is potentially better than that of any single device. In practice, the measurements from many inexpensive MEMS devices on a single chip can be analytically combined in a real-time processing filter to raise their performance from low-grade to high grade, and command a much wider range of applications.

[0009] The basic prior art approach of Allan was to study the behavior of the three sports watches relative to coordinated universal time (UTC) time, and relative to each other. This formed a basis for combining them into an improved time estimate. The present approach is different in at least two important ways. First, gyro rate measurements are combined rather than time measurements. Second, a minimum-variance estimation approach and Kalman filtering are used to provide a systematic theory and implementation for combining devices. The Kalman filtering approach has the advantage of being a systematic method to ensure a minimum variance rate estimate.

[0010] The invention utilizes the collective behavior of multiple MEMS devices to offset deficiencies occurring at the single-device level, thereby improving overall performance to the point where they can be considered as a competitive technology for future NASA/JPL missions, aircraft, land vehicles and underwater craft.

[0011] The main theoretical result is a complete optimal minimum variance method for combining and filtering measurements from multiple gyro devices. An optimal filter combines and smoothes multiple gyro measurements to produce a minimum-variance rate estimate. This rate estimate is denoted as a “virtual gyro” since it is indistinguishable from an actual rate measurement taken from a single physical gyro having improved performance.

[0012] The gain matrix of the minimum-variance filter is an important quantity because it prescribes how the measurements are weighted and smoothed with respect to each other. According to Kalman theory, the gain matrix of the optimal filter is uniquely specified by the filter covariance which propagates according to a matrix Riccati differential equation (RDE). An important contribution of this research effort is the complete and exact analytical solution to the virtual gyro Riccati differential equation. This takes a large step toward understanding virtual gyros and being able to synthesize them in practice. The analytic solution reveals a wealth of theoretical properties, and leads, for example, to a simple parallel resistor analogy for understanding how the component gyro drifts add to form the overall virtual gyro drift. The analytic solution also leads to an important discrete-time implementation of the minimum variance filter.

[0013] For example, an interesting case study is given where four gyros, each with 8.6 degrees per hour drift, are combined to give a virtual gyro drift of 0.062 degrees per hour. This indicates a factor of over 100 times improvement. While performance is generally problem specific, results indicate that substantial improvement is possible when there are favorable correlations between the component devices.

[0014] More specifically the invention is defined as a method of processing signals from multiple gyros each having a rate, in order to derive a higher accuracy output signal. The method of processing signals comprising the steps of measuring the rate of multiple gyro signals, and combining the measured rates using a recursive processing filter whose properties are optimized using the statistical correlations and noise characteristics of the multiple rate sensing devices. A specific method for recursive processing is based on a steady state Kalman filter to generate a minimum variance rate estimate. Consider N single-axis gyros, each of which provides a noisy measurement yi of the true angular rate ω, i.e.,

y i =ω+b i +n 1i ; i=1, . . . , N

[0015] where bi and n1i denote noises associated with the ith device. Specifically, the term n1i is a delta-correlated noise process on the ith device, called the angle-random walk (ARW). The term bi is a slowly varying random quantity called the gyro bias. The gyro bias for the ith device is modeled by,

{dot over (b)}i=n2i

[0016] where the term n2i is a delta-correlated noise process called the rate-random walk (RRW). Let 1 denote an N vector completely containing ones, i.e., 1$1\stackrel{\Delta }{=}\left[\begin{array}{c}1\\ ⋮\\ 1\end{array}\right]$

[0017] With this notation, it is convenient to “stack” the above relations and write the collection of measured signals in the following vectorized form,

y=b+ 1ω+v

{dot over (b)}=wb

[0018] where y=[y1,y2, . . . , yN]T is an N vector comprised of the measured signals (i.e., the virtual gyro input signals) from N gyros; b=[b1, . . . , bN]T is the vector of gyro biases; ω is the true rate (common to all gyro devices); v=[n11, n12, . . . , n1N]T is a vector of ARW component noises; and wb=[n21, n22, . . . , n2N]T is a vector of RRW component noises. The true rate ω is modeled as a random walk of the form,

{dot over (ω)}=wω

[0019] where wω is a delta-correlated noise process.

[0020] An N+1 dimensional state vector will be defined as, 2$x=\left[\begin{array}{c}b\\ \omega \end{array}\right]$

[0021] so that the measurement vector y can be written equivalently as,

y=Hx+v

[0022] where H is the N×(N+1) matrix given by,

H=[I1]

[0023] and where I denotes the N×N identity matrix. It is also convenient to define a N+1 selection vector eN+1 that contains all zeros except for a one in the last place, i.e., 3$e_{N+1}=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 1\end{array}\right]$

[0024] The selection vector is useful for extracting the rate signal from the state vector as, 4$\omega =e_{N+1^x}^T$

[0025] Two important quantities are the covariance matrices R and Q which are defined as follows,

[0026] RCov[v] is an N×N covariance matrix of the ARW noise. The matrix R can have off-diagonal terms to account for correlations between the component ARW noises of different gyros. 5$Q\stackrel{\Delta }{=}\left[\begin{array}{cc}{Q}_b& 0\\ 0& q_{\omega }\end{array}\right]$

[0027] is an augmented RRW matrix, where QbCov[wb] is the N×N covariance matrix of the RRW noise wb, and qω=Cov[wω] is the scalar covariance of the delta-correlated noise wω used in the model of the true rate ω. The matrix Qb can have off-diagonal terms to account for correlations between the component RRW noises of different gyros.

[0028] Several key matrix quantities are computed directly from the ARW matrix R and the augmented RRW matrix Q. Specifically, these are the (N+1)×(N+1) matrices A, D and S, the (N+1)×N matrix {overscore (A)}, and the N×N matrix {overscore (D)}. Details on how these matrices are computed are treated in detail in Bayard, and will be summarized in the detailed invention description section. However, an important new innovation is that the matrices A,{overscore (A)},D,{overscore (D)},S can be computed analytically without having to numerically integrate a nonlinear matrix Riccati equation. This represents a sharp break with standard Kalman filtering approaches, yet is critical for realizing the minimum-variance filter in practice, since the Riccati equation diverges with time making it impossible to otherwise implement a steady-state filter using finite-precision arithmetic.

[0029] Once the matrices A,{overscore (A)},D,{overscore (D)},S are available, the gain matrix K∞ for the continuous-time minimum-variance filter and the gain matrix {overscore (K)} for the discrete-time minimum-variance filter can be computed as follows, 6$K_{\infty }=Q^{1/2}\stackrel{\_}{A}{\stackrel{\_}{D}}^{-1/2}{\stackrel{\_}{A}}^T{Q}^{1/2}{H}^T{R}^{-1}$$\stackrel{\_}{K}=Q^{1/2}\stackrel{\_}{A}{\stackrel{\_}{D}}^{-1/2}\left(I-e^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}\right){\stackrel{\_}{A}}^T{Q}^{1/2}{H}^T{R}^{-1}$

[0030] Here, the notation X−1, X1/2, XT, eX denotes the matrix inverse, matrix square-root, matrix transpose, and matrix exponential, respectively, for an arbitrary square positive definite or semi-definite matrix X.

[0031] The gain matrices K∞ and {overscore (K)} are used in conjunction with the A,{overscore (A)},D,{overscore (D)},S matrices to give the continuous-time and discrete-time filter implementations summarized below.

[0032] Continuous-Time (i.e., Analog) Implementation of Minimum-Variance Filter Continuously-sampled analog measurements from multiple gyros, denoted as y(t), are processed in a minimum-variance filter defined by the following system of differential equations, 7$\stackrel{.}{\hat{\xi }}=\left[\begin{array}{cc}-{\stackrel{\_}{D}}^{\frac{1}{2}}& 0\\ 0& 0\end{array}\right]\hat{\xi }+S^{-1}{K}_{\infty }y$$\hat{\omega }=e_{N+1}^{T}S\hat{\xi }$

[0033] where K∞ is the minimum-variance gain matrix, {circumflex over (ξ)} is the continuous-time filter state vector of dimension N+1, {circumflex over (ω)} is the optimized minimum-variance rate estimate (i.e., the resulting virtual gyro output signal).

[0034] Comments

[0035] The filter state {circumflex over (ξ)} is a minimum-variance estimate of ξ which is defined by the transformation ξ=S−1x of the true state x. This transformation is important because it explicitly isolates the pure integral action so that practical measures can be taken to ensure that it is not destabilized by numerical round-off errors in the implementation. Two such practical methods are given below denoted as (i) and (ii), without which the minimum-variance filter would not be realizable with finite-precision arithmetic.

[0036] (i) The last row of the matrix product S−1 K∞ will be exactly zero, according to the new analytical solution. Accordingly, the implementation requires that the last row of S−1 K∞ (after it is computed arithmetically) be explicitly set to zero. This avoids destabilizing effects due to imperfect finite-arithmetic errors in the computation.

[0037] (ii) The isolated “0” in the lower right-hand corner of the matrix 8$\&AutoLeftMatch;\left[\begin{array}{cc}-{\stackrel{\_}{D}}^{\frac{1}{2}}& 0\\ 0& 0\end{array}\right]$

[0038] in the filter implementation should be set to zero using the full numerical precision available on the machine or apparatus implementing the filter, to ensure (using finite-precision arithmetic) that for any finite-precision number a the result of multiplication is a·0=0.

[0039] Discrete-Time (i.e., Digital) Implementation of Minimum Variance Filter

[0040] Time-sampled measurements from multiple gyros, denoted as y(kT) for the kth time instant, are processed in a minimum-variance filter given by the following recursive equations, 9$\hat{\xi }\left(\left(k+1\right)T\right)=\left[\begin{array}{cc}{e}^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}& 0\\ 0& 1\end{array}\right]\hat{\xi }\left(kT\right)+S^{-1}\stackrel{\_}{K}y\left(kT\right)$$\hat{\omega }\left(kT\right)=e_{N+1}^{T}S\hat{\xi }\left(kT\right)$

[0041] where {overscore (K)} represents the minimum variance gain matrix; T is the sampling period, {circumflex over (ξ)} is the discrete-time filter state vector of dimension N+1, {circumflex over (ω)}(kT) is the optimized minimum-variance rate estimate (i.e., the resulting virtual gyro output signal) at time instant kT.

[0042] Comments

[0043] The discrete-time filter state {circumflex over (ξ)}(kT) is a minimum-variance estimate of ξ(kT) which is defined by the transformation ξ(kT)=S−1x(kT) of the true state x. This transformation is important because it explicitly isolates the pure integral action so that practical measures can be taken to ensure that it is not destabilized by numerical round-off errors in the implementation. Two such practical methods are denoted as (i) and (ii) below, without which the minimum-variance filter would not be realizable with finite-precision arithmetic.

[0044] (i) The new theory dictates that the last row of the matrix product S−1 {overscore (K)} will be exactly zero. Accordingly, the implementation requires that the last row of S−1 {overscore (K)} (after it is computed arithmetically) be explicitly set to zero. This avoids destabilizing effects due to imperfect finite-arithmetic errors in the computation.

[0045] (ii) The isolated “1” seen in the lower right-hand corner of the matrix 10$\&AutoLeftMatch;\left[\begin{array}{cc}{e}^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}& 0\\ 0& 1\end{array}\right]$

[0046] should be set to unity using the full numerical precision available on the machine implementing the filter, to ensure (using finite-precision arithmetic) that for any finite-precision number a the result of multiplication is a·1=a.

[0047] The step of combining the measured rates of the multiple gyros using a steady state Kalman filter to generate a minimum variance rate estimate can also be said to comprise the steps of: sampling N analog signals from the N gyros over a period of T seconds; digitizing the sampled N signals into a digital data point for each sampling period T for each gyro; assembling the N digitized signals from the N gyros into the N-vector y(kT); processing y(kT) by multiplying it with the matrix {overscore (K)} to generate the product signal {overscore (K)}y(kT); processing the {overscore (K)}y(kT) signal by multiplying it with a matrix S−1 to generate the product signal S−1 {overscore (K)}y(kT); inputting the S−1 {overscore (K)}y(kT) product signal during a kth sampling period into a first input of a summing node having two inputs; outputting the sum of the two inputs of the summing node to a delay operator; delaying the contents of the delay operator by one sampling period, T; feeding back the delayed output of the delay operator to a smoothing matrix operator modeled by 11$\left[\begin{array}{cc}{e}^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}& 0\\ 0& 1\end{array}\right]$

[0048] to generate a feedback signal; inputting the feedback signal to a second one of the inputs to the summing node such that it is added to the S−1 {overscore (K)}y(kT) product signal generated in the k+1 sampling period; processing the contents of the delay operator for the kth sampling period by multiplying it with a matrix S operator to generate a vector output signal; and extracting the last matrix element of the vector output signal as the best estimate of the virtual gyro rate ŵ(kT) for the kth sampling point.

[0049] The invention is also to be understood as comprising an apparatus for performing each of the foregoing steps. The means for performing the steps may be comprised of appropriately designed logic circuits, firmware controlled digital processors and/or software controlled computers in any combination desired among the various steps and means.

[0050] While the invention has been specifically described in terms of gyros, the invention is fully and equally applicable to any plurality of devices which generate a rate or other drifting signal subject to noise. Such devices may include clocks, accelerometers, or any other system of multiple sensors intended to measure the same physical quantity over time. The correlation among the plurality of devices is entirely general and may assume any type of correlated pattern now known or later devised.

[0051] While the apparatus and method has or will be described for the sake of grammatical fluidity with functional explanations, it is to be expressly understood that the claims, unless expressly formulated under 35 USC 112, are not to be construed as necessarily limited in any way by the construction of “means” or “steps” limitations, but are to be accorded the full scope of the meaning and equivalents of the definition provided by the claims under the judicial doctrine of equivalents, and in the case where the claims are expressly formulated under 35 USC 112 are to be accorded full statutory equivalents under 35 USC 112. The invention can be better visualized by turning now to the following drawings wherein like elements are referenced by like numerals.

BRIEF DESCRIPTION OF THE DRAWINGS

[0052] FIG. 1 is a block diagram of a mathematical model by which the Kalman filtering method of the invention is illustrated.

[0053] FIG. 2 is a block diagram of a continuous time model of the Kalman filtering method of the invention.

[0054] FIG. 3 is a block diagram of a discrete time model of the Kalman filtering method of the invention.

[0055] FIG. 4 is a graph of the normalized drift of a four component virtual gyro verses different correlation factor values.

[0056] FIG. 5 is a graph of the drifting rate measurements from four uncorrelated component gyros verses time.

[0057] FIG. 6 is a graph of the optimal minimum variance rate of a virtual gyro using the component gyros of FIG. 5.

[0058] FIG. 7 is a graph of the drifting rate measurements from four favorably correlated component gyros verses time.

[0059] FIG. 8 is a graph of the optimal minimum variance rate of a virtual gyro using the component gyros of FIG. 7.

[0060] FIG. 9 is a graph of the drifting rate measurements from four very favorably correlated component gyros verses time.

[0061] FIG. 10 is a graph of the optimal minimum variance rate of a virtual gyro using the component gyros of FIG. 9.

[0062] The invention and its various embodiments can now be better understood by turning to the following detailed description of the preferred embodiments which are presented as illustrated examples of the invention defined in the claims. It is expressly understood that the invention as defined by the claims may be broader than the illustrated embodiments described below.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0063] Gyros are described in the illustrated embodiment rather than accelerometers, since gyros benefit most from new methods for reducing drift. Consider first the use of an optimal minimum-variance filter. The main result is a complete optimal minimum variance method for combining and filtering multiple gyro measurements. A mathematical model for the methodology is shown in FIG. 1. An angle random walk (ARW) noise source 12 provides input signals through a mixing matrix 16 to a plurality of summing points 14(1), 14(2), . . . 14(N) for N gyros. A rate random walk (RRW) noise source 20 similarly provides input signals through a mixing matrix 18 to a plurality of summing points 14(1), 14(2), . . . 14(N) for the N gyros. The signals from matrix 18 are integrated in integrators 22(1), 22(2) . . . 22(N) and 24, the last of which corresponds to a true or perfect drift-free ideal gyro providing the true rate. This truth model is a mathematical construct and physically does not exist. The true and ARW and RRW noise signals are summed at points 14(1), 14(2), . . . 14 (N) and output as multiple gyro measurements going into the Kalman filter 10. The output of filter 10 is the desired minimum variance rate estimate for the N gyros. For illustration only (not used by the filter) the output of filter 10 is combined at summing point 26 with the negative true rate to define the rate error, which is the mathematical quantity of interest that is being minimized by this overall approach in a minimum-variance sense.

[0064] Here, an optimal filter 10 is provided to combine and smooth multiple gyro measurements (all oriented along the same axis) to produce a minimum-variance rate estimate (i.e., the mean-square value of the quantity “Rate Error” in FIG. 1 at output 28 is minimized by choice of this filter 10). Even though this rate estimate is analytically or mathematically derived from N gyros measurements as shown in FIG. 1, it is indistinguishable from a single rate measurement which might be taken from a single physical gyro. Consequently the resulting optimal minimum-variance rate estimate is denoted as a “virtual gyro” measurement.

[0065] Each gyro is considered to provide a signal yi=ω+bi+n1i for i=1, . . . , N, where yi is the measured angular rate of the ith gyro, ω is the true rate, bi is a slowly varying random quantity of the ith gyro denoted as the gyro bias, n1i is a delta-correlated noise process of the ith gyro denoted as the angle-random walk ARW noise, and the vector of ARW noises is denoted as v=[n11, n12, . . . , n1N]T. The measurement equation for the N gyro system then can be written in vector form as y=b+1ω+v, where each of the quantities are column vectors. The N gyros, however, are not independent, but have a drift or noise which is correlated because of the physical correlations of their manufacture and design. Hence, we can assign a covariance matrix R for ARW and Q for RRW representing the degree of noise correlation between each of the gyros relative to themselves and relative to each of the other gyros. Greater detail concerning these correlation matrices are given in the attached appendix.

[0066] The gain matrix K contained inside the minimum-variance filter 10 is an important quantity because it prescribes how the measurements are weighted with respect to each other, and completely characterizes how the dynamics of the filter 10 are assigned to optimally smooth the data. According to Kalman theory, the gain matrix K of the optimal filter 10 is a function of the (N+1)×(N+1) state covariance matrix P(t) which propagates according to a matrix Riccati differential equation (RDE) (See Bucy et.al., “Filtering for Stochastic Processes with Application to Guidance,” Wiley, New York 1968). The gain matrix K and optimal filter associated with the virtual gyro is characterized by the set of matrix equations below: 12$\begin{array}{cc}\stackrel{\stackrel{.}{^}}{x}=K\left(t\right)\left(y-H\hat{x}\right)\stackrel{.}{P}=Q-P{H}^T{R}^{-1}HPK\left(t\right)=P\left(t\right)H^T{R}^{-1}& \left(1\right)\end{array}$

[0067] where {dot over (P)} is the time derivative of the state covariance matrix P(t); Q is an (N+1)×(N+1) augmented RRW covariance matrix; H is an N×(N+1) matrix in the measurement equation y=Hx+v, where the ARW noises have been stacked into the vector v=[n11, n12, . . . , n1N]T; y is a vector representing the measured signals from N gyros, x is the state vector x=[bT,ω]T; {circumflex over (x)} is a minimum-variance estimate of the state vector x; R−1 is the inverse of the N×N ARW covariance matrix R; and K(t) is the time dependent Kalman gain.

[0068] In general, the Riccati differential equation is a nonlinear matrix differential equation, and is very difficult to solve analytically. Because of this, solutions used in practice are based on integrating the differential equations numerically, or by calculating a steady-state solution using iterative numerical methods.

[0069] Unfortunately, there are two complications associated with the virtual gyro Riccati differential equation which prevent its solution using these standard approaches. The first problem is that a steady-state solution to the Riccati differential equation does not exist due to lack of “complete observability” as described in greater detail in the appendix (See also Bucy). Physically, this important property does not hold because the optimal estimate, however good, always has some non-zero component of drift.

[0070] The second problem, is that the presence of drifts causes the covariance P(t) to increase over time so that the numerical Riccati differential equation solution is unstable and becomes unbounded in the limit of time t going to infinity. Because of this divergence with time, solving for K(t) from equation (1) using finite-precision arithmetic tends to result in an unbounded or poorly behaved solution which has no physical meaning. This is because with finite-precision arithmetic, the increase in P cannot be perfectly cancelled by post-multiplication by HT and hence the product typically becomes large as time t goes to infinity. In contrast, the analytical solution of the Riccati differential equation for K(t) provides a method whereby a finite and meaningful solution can be obtained.

[0071] As shown in the attached appendix it is extremely fortunate that certain covariance terms vanish after post-multiplying by HTR−1. This implies that the Kalman gain matrix approaches a steady-state value with time, even though the covariance from the RDE solution diverges. The resulting constant gain matrix K∞ (associated with the continuous-time filter implementation), is given by any of the following expressions, 13$\begin{array}{cc}{K}_{\infty }=Q^{\frac{1}{2}}{L}^{†/2}{Q}^{\frac{1}{2}}{H}^T{R}^{-1}=Q^{\frac{1}{2}}A{D}^{†/2}{Q}^{\frac{1}{2}}{H}^T{R}^{-1}=Q^{\frac{1}{2}}\stackrel{\_}{A}{\stackrel{\_}{D}}^{-\frac{1}{2}}{\stackrel{\_}{A}}^T{Q}^{\frac{1}{2}}{H}^T{R}^{-1}& \left(2\right)\end{array}$

[0072] where L=Q1/2HTR−1HQ1/2, the symbol “†” denotes the Moore-Penrose pseudoinverse (which always exists), and the matrices A,{overscore (A)},D,{overscore (D)} are defined by an eigenvalue decomposition of the symmetric matrix L (discussed below, and in the appendix).

[0073] It is emphasized that these formulas allow the optimal filter to be determined completely algebraically. This avoids the need for numerical integration of the RDE, and is therefore not subject to rounding errors incurred from inexact cancellation of any unbounded term when forming the Kalman gain. Starting from equation (2) it can be shown that in the digital case, the equivalent constant filter gain K can be directly calculated using 14$\begin{array}{cc}\stackrel{\_}{K}=S\left[\begin{array}{cc}{\stackrel{\_}{D}}^{-\frac{1}{2}}\left(I-e^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}\right)& 0\\ 0& 0\end{array}\right]S^T{H}^T{R}^{-1}=Q^{1/2}\stackrel{\_}{A}{\stackrel{\_}{D}}^{-1/2}\left(I-e^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}\right){\stackrel{\_}{A}}^T{Q}^{1/2}{H}^T{R}^{-1}& \left(3\right)\end{array}$

[0074] where S=Q1/2A. Thus a major contribution of this research effort has been the complete and exact analytical solution to the virtual gyro Riccati matrix differential equation. This result is summarized at length and in detail in Bayard et.al. “Foundations of Virtual Gyro Synthesis” JPL D-21656 Jan. 21, 2002, which is attached and incorporated herein as an appendix, which was filed as part of the original application and is now available to the public in the file wrapper of this application, but which has been deleted from the printed specification for the sake of brevity.

[0075] We consider this result a major accomplishment in the field, and a large step forward in understanding virtual gyros and being able to synthesize them in practice. The analytic solution reveals a wealth of theoretical properties, and enables for the first time the consideration of practical implementations.

[0076] Many important properties of virtual gyros become transparent as a consequence of knowing the analytic solution to the underlying Riccati differential equation. These properties have been rigorously derived and documented in Bayard. The main highlights and their significance are briefly discussed below.

[0077] Even though the Riccati differential equation grows in an unbounded fashion as a function of time, the Kalman gain K(t) can be shown to approach a steady-state matrix. This result is both remarkable and fortuitous since it significantly simplifies implementations and, for the first time, enables practical schemes for realizing the optimal filter using a constant minimum-variance gain matrix K∞.

[0078] The analytic solution to the Riccati differential equation allows the development of a complete statistical theory which characterizes virtual gyro drift, and provides theoretical limits of improvement due to using an ensemble of correlated sensing devices as a single virtual sensor.

[0079] The transfer function of the minimum-variance filter is denoted by K(s) (in the Laplace variable s). K(s) has been analyzed in detail and its properties revealed as shown in the appendix. The transfer function K(s) indicates the presence of a marginally unstable pole, which would otherwise make computational implementation impossible if not properly understood and compensated. The analytic-based implementation approach provides a method to compensate for the unstable pole. In addition, a simple algebraic method has been developed for computing the optimal gain matrix K∞. With this result, the finite-precision arithmetic issues associated with the forming the gain matrix from the Riccati differential equation can be avoided completely in the implementation.

[0080] The notion of statistical common mode rejection has been characterized mathematically. For statistically uncorrelated gyros (i.e., where the Q and R matrices are diagonal), it has been shown in Bayard that component drifts add like parallel resistors in terms of their variances (e.g., in units of rad2/sec3). For example, the drift of a virtual gyro composed of N independent and identical devices will be 1/N of the individual device variance e.g., like N resistors R in parallel having overall resistance R/N. Note this corresponds to 1/{square root}{square root over (N)} of the individual device drift in units of rad/sec3/2. Along the same lines, any single device with a large drift will not count much in the overall performance (like a large resistor among a parallel connection of resistors, which can be removed from consideration as a practical matter). Conversely, any device having drift significantly less than any of the other devices will dominate the performance of the virtual gyro (like a small resistance among a parallel connection of resistors).

[0081] Potential improvement in drift is much more impressive when the noise between devices is correlated (i.e., where the Q and R matrices have off-diagonal elements). Specifically, a key result in Bayard proves that the theoretically minimum drift attainable from combining N gyro devices having Rate Random Walk matrix Qb is given by, 15$D_{\mathrm{vg}}=\frac{1}{\mathrm{Sum}\mathrm{over}\mathrm{elements}\mathrm{of}{Q}_b^{-1}}$

[0082] The units of Dvg are rad2/sec3 which can be transformed to an equivalent one-sigma random walk growth in units of deg/hr over a one hour period. Drift values reported in the remainder of this document in units of deg/hr will have this specific interpretation. It is important to note that this formula gives the theoretically lowest value of gyro drift that is attainable by combining measurements from multiple gyro devices. In this light, it is emphasized that a real-time recursive processing filter which is theoretically capable of attaining this minimum drift bound is the subject of the present invention.

[0083] For correlated devices, interpretations given in Bayard indicate that the “principal component” drifts divided by their common-mode device fraction (a measure of how coherently the noise sources influence the multiple devices) add like parallel resistors. Roughly speaking, this says that the component drifts count to the extent that they influence multiple devices coherently (e.g., with common sign), and can be removed to the extent that they influence multiple devices incoherently (e.g., with randomized signs). This generalizes the intuitive notion that a noise which is additive to two devices but with opposite sign, can be eliminated by simple averaging. Noise eliminated by such favorable correlations can improve performance of the virtual gyro far beyond the 1/{square root}{square root over (N)} factor attainable using uncorrelated devices.

[0084] The main significance of the above results, is that virtual gyros can now be developed and implemented using N physical gyros. The present theory optimally combines individual sensing devices, and offers a method to improve overall performance characteristics beyond what is possible in any single device. A rough estimate, assuming that the correlation between devices in close spatial proximity is favorable, is that between one or two orders of magnitude improvement in drift is possible.

[0085] Consider now various practical virtual-gyro implementations. The procedure to compute the gain matrix K∞ for the continuous-time minimum-varance filter implementation, and the gain matrix {overscore (K)} for the discrete-time minimum-variance filter is now given. The covariance matrices Q and R are assumed known from either engineering and/or testing of the gyros, because how the drift in one gyro is related to another can be determined empirically or by first principles. According to the invention, the filter which provides the best estimated value of the gyro rate by combining the N component gyro measurements is achieved when the minimum-variance filter gain K∞ (for the continuous-time problem) and {overscore (K)} (for the discrete-time problem) is known. The steps for determining these gains from the Q and R matrices is given below.

[0086] Step (1) Form the (N+1)×(N+1) symmetric matrix L=Q1/2HTR−1HQ1/2

[0087] Step (2) Perform the eigenvalue decomposition of L ordered in such a manner that its zero eigenvalue is last,

L=ADA T

[0088]  where A is an orthogonal matrix, and D is a diagonal matrix whose last diagonal element is zero. As a practical matter the smallest eigenvalue will never be exactly zero when using finite-precision arithmetic, in which case the smallest eigenvalue should be used to represent zero for the purposes of the eigenvalue ordering.

[0089] Step (3) Extract the N×N matrix {overscore (D)} and (N+1)×N matrix {overscore (A)} from D and A, respectively, as follows, 16$D=\left[\begin{array}{cc}\stackrel{\_}{D}& 0\\ 0& 0\end{array}\right]$$A=\left[\begin{array}{ccc}\stackrel{\_}{A}& ⋮& v\end{array}\right]$

[0090] Step (4) Form the auxiliary (N+1)×(N+1) matrix S, and the (N+1)×1 selection vector eN+1 as follows, 17$S=Q^{1/2}A$$e_{N+1}^T=\left[0,\dots ,0,1\right]$

[0091] Step (5) Form the Kalman gain for a continuous-time implementation as,

K ∞ =Q 1/2 {overscore (A)}{overscore (D)} −1/2 {overscore (A)} T Q 1/2 H T R −1

[0092] Step (6) Form the Kalman gain for a discrete-time (digital) implementation as, 18$\stackrel{\_}{K}=S\left[\begin{array}{cc}{\stackrel{\_}{D}}^{-1/2}\left(I-e^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}\right)& 0\\ 0& 0\end{array}\right]S^T{H}^T{R}^{-1}=Q^{1/2}\stackrel{\_}{A}{\stackrel{\_}{D}}^{-1/2}\left(I-e^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}\right){\stackrel{\_}{A}}^T{Q}^{1/2}{H}^T{R}^{-1}$

[0093] Continuous-Time Implementation

[0094] The continuous-time filter implementation will be considered first. Using the Kalman gain, the virtual gyro rate estimate can be propagated in state-space form as, 19$\stackrel{\stackrel{.}{^}}{x}=K_{\infty }(y-\hat{y})=-K_{\infty }H\hat{x}+K_{\infty }y$$\hat{\omega }=e_{N+1}^T\hat{x}$

[0095] where the variables with the “hat” are the statistical minimum variance estimates. The transfer function from the measurements to the virtual gyro estimate can be written in the Laplace transform domain as, 20$\hat{\omega }\left(s\right)={e_{N+1}^T\left(sI+K_{\infty }H\right)}^{-1}{K}_{\infty }y\left(s\right)$

[0096] where for this equation the quantities {circumflex over (ω)} and y are to be interpreted as Laplace-transformed variables. Unfortunately, this filter has a pole exactly at the origin. This implementation will not be reliable due to numerical rounding errors that can cause the pole to become unstable. A better approach is developed below. A main emphasis in the appendix has been to study the optimal filter matrix K∞H closely, since it completely characterizes the poles of the Kalman filter.

[0097] We turn now to the filter implementation. After some mathematical manipulation as set forth in detail in the appendix the continuous time optimal filter 10 of FIG. 2 can be realized. The optimal gain K∞ is a function of the drift rate signal model parameter qω, which can be interpreted as a tuning parameter in this filter design. One might think of making qω large to minimize any dependence of the filter on the a-priori smoothness of the true rate. While theoretically possible, the poles of the optimal continuous-time filter approach negative infinity and become difficult to implement as qω becomes large. In contrast, the poles of the digital filter considered below move towards the origin in the limit as qω becomes large, and approach an easily implemented deadbeat condition. FIG. 2 is a diagram that illustrates the virtual gyro implementation using a continuous-time minimum variance filter. The last row of the input matrix S−1 K∞ (corresponding to the pure integrator state) should be set exactly to zero. This is dictated by the theory but must be explicitly enforced because it otherwise may not be so due to finite-precision arithmetic considerations. While the continuous-time filter can be implemented using analog circuitry, it is becoming more convenient to utilize computers and digital circuitry so that a digital implementation is considered next.

[0098] Digital (Discrete-Time) Implementation

[0099] A digital implementation of the virtual gyro is considered in the case where individual gyro devices are sampled in time with sampling period T. FIG. 3 is a diagrammatic virtual gyro implementation using a discrete-time minimum variance filter. The matrix product S−1 {overscore (K)} can be formed directly, or by using the alternative expression derived in the appendix, 21$\begin{array}{cc}{S}^{-1}\stackrel{\_}{K}=S\left[\begin{array}{cc}{\stackrel{\_}{D}}^{-\frac{1}{2}}\left(I-e^{-{\stackrel{\_}{D}}^{\frac{1}{2}}T}\right)& 0\\ 0& 0\end{array}\right]S^T{H}^T{R}^{-1}& \left(4\right)\end{array}$

[0100] The last row of the input matrix S−1 {overscore (K)} (corresponding to the pure integrator state) should be set exactly to zero.

[0101] FIG. 3 is a clocked digital system. The analog signal from each gyro is sampled and digitized into a digital number or data point for each sampling period T. The N signals from the N gyros are assembled at buffer 30 into an N-vector which is then fed into a matrix multiplier 32 of first the {overscore (K)} operator 32a and then the matrix S−1 operator 32b, for which the combined operation is given by the above matrix equation (4). The matrix {overscore (K)} is of dimension (N+1)×N and the matrix S−1 operation 32b is of dimension (N+1)×(N+1) so that the input to the S−1 {overscore (K)} operator 32 is an N vector, while the output is an N+1 vector. The matrix product result is clocked into summing node or adder 34. The output of summing node 34 is clocked into the operator z−1 I (where I is the identity matrix, and z−1 is the one-sample delay) which simply applies a one clock delay, which is then fed back through a smoothing matrix 22$\left[\begin{array}{cc}{e}^{-{\stackrel{\_}{D}}^{\frac{1}{2}}}& 0\\ 0& 1\end{array}\right]$

[0102] operator 38 to summing node 34. At this point the modified value of the prior data point is added to the next data point. Meanwhile the previous addition result is forwarded and multiplied by the matrix S operator 40. The output of the matrix S operator 40 is multiplied by the 23$e_{N+1}^T$

[0103] matrix operator 42. The 24$e_{N+1}^T$

[0104] matrix is all zeroes except for the last row and column, which is 1. This 25$e_{N+1}^T$

[0105] matrix operator 42 simply picks out the last matrix element of its input which is the desired answer, i.e., the best estimate of the virtual gyro rate {circumflex over (ω)}(kT) for the kth sampling point. Elements 30-42 have been described as buffers or operators, but they are to be understood to comprise any type of signal processing circuit including custom logic circuits, firmware controlled digital processors, or software controlled computers.

[0106] A critical observation (analogous to the continuous-time case) is that the last row of the input matrix S−1 {overscore (K)} is zero, so that the integrator state is unforced and will remain constant for all time. In practice, the designer should set the last row of S−1 {overscore (K)} exactly to zero, since numerical rounding errors may otherwise cause it to have small nonzero entries. This approach completely resolves any potential stability problems arising from implementing the optimal filter.

EXAMPLES

[0107] Some simulation examples will demonstrate the virtual gyro implementation using the optimal digital filter depicted in FIG. 3. All simulations use N=4 gyros, with a sampling rate of T=1 sec, and a total running time of 400,000 secs (approximately 111 hours). The three examples correspond to the choices of the cross correlation factor ρ being set equal to 0; −0.32890; and −0.33329 where it is assumed that all gyros have constant cross-correlation factor ρ. The normalized virtual gyro drift can be calculated by setting N=4 in (5.10) to yield a theoretical expression that predicts that a perfect driftless virtual gyro is obtained in the limit as the correlation factor p approaches −⅓. The choices of p are intentionally chosen to successively approach this limiting value. They are plotted in FIG. 4 with a “*”, which shows the normalized virtual gyro drift versus −ρ and the expected drift reduction for each case. Specifically, the expected improvement for each of the choices of ρ can be calculated based on the analytical solution, and then compared to a computer simulation. The case ρ=0 corresponds to uncorrelated gyro devices for which one expects an improvement of {square root}{square root over (N)}={square root}{square root over (4)}=2. The choice ρ=−0.32890 corresponds to favorably correlated gyros with an improvement of 17.3 times, and the choice ρ=−0.33329 corresponds to very favorably correlated gyros with an improvement of 173 times.

[0108] Consider first the uncorrelated gyro case. The first example considers the simulated choices N=4, qb=6.5290e−013 (rad2/sec3), and ρ=0. Since ρ=0, the gyros are uncorrelated. Here, each gyro device by itself has a drift of 10 (deg/hr). Since the four gyros are uncorrelated by design, one would expect a drift reduction factor of {square root}{square root over (N)}={square root}{square root over (4)}=2 or equivalently 10/2=5 (deg/hr) in the final virtual gyro design. The simulated time responses of the four gyro devices are shown in FIG. 5. The resulting virtual gyro response is shown in FIG. 6. The drift is reduced from 10 (deg/hr) for the individual gyros to a value of 4.9301 (deg/hr) for the combined gyros. This agrees well with the theoretically predicted value of 5 (deg/hr).

[0109] In the case of favorably correlated gyros the choices N=4, qb=4.9131e−013 (rad2/sec3), and ρ=−0.32890. The gyro devices are intentionally correlated to give an expected theoretical factor of 17.3 magnitude reduction in drift when combined. Here, each gyro device by itself has a drift of 8.6747 (deg/hr). The theoretical virtual gyro drift can be calculated to give, 1.6350e−015 (rad2/sec3). This corresponds to a drift of 0.5 (deg/hr) in one hour, which is a 17.3 times magnitude improvement over the 8.6747 (deg/hr) component gyro drifts.

[0110] The time responses of the four gyro devices are simulated and shown in FIG. 7. The resulting virtual gyro response is shown in FIG. 8. The drift is reduced to a value of 0.51883 (deg/hr) by combining the devices. This is over an order of magnitude improvement relative to the individual gyro variances, and agrees well with the 0.5 (deg/hr) expected theoretical drift rate.

[0111] Consider finally a case of very favorably correlated gyros. Here we choose N=4, qb=4.8969e−013 (rad2/sec3), and ρ=−0.33329 in (5.9). As discussed earlier, the gyro devices are intentionally correlated to give a factor of 173 times magnitude reduction in drift when combined. The RRW and ARW matrices are given specially as, 26$Q_b=\left[\begin{array}{cccc}4.8969e-13& -1.632e-13& -1.632e-13& -1.632e-13\\ -1.632e-13& 4.8969e-13& -1.632e-13& -1.632e-13\\ -1.632e-13& -1.632e-13& 4.8969e-13& -1.632e-13\\ -1.632e-13& -1.632e-13& -1.632e-13& 4.8969e-12\end{array}\right]\left(\frac{{\mathrm{rad}}^2}{{\sec }^3}\right)$$R=\left[\begin{array}{cccc}7.6366e-10& 0& 0& 0\\ 0& 7.6366e-10& 0& 0\\ 0& 0& 7.6366e-10& 0\\ 0& 0& 0& 7.6366e-10\end{array}\right]\left(\frac{{\mathrm{rad}}^2}{{\sec }^5}\right)$

[0112] Here, each gyro device by itself has a drift of 8.6604 (deg/hr). The theoretical virtual gyro drift can be calculated by using N=4; qb=4.8969e−013; ρ=−0.33329 to give 1.6e−017 (rad2/sec3). This corresponds to a expected theoretical drift of 0.05 (deg/hr) in one hour, which is a 173 times magnitude improvement over the 8.6604 (deg/hr) component gyro drifts.

[0113] The time responses of the four gyro devices are simulated and shown in FIG. 9. A virtual gyro is designed using the digital observer approach shown in FIG. 3, with the choice of qω=4.8414e−009 (rad2/sec3) for rate signal modeling. The resulting virtual gyro response is shown in FIG. 10. The drift is reduced by combining the devices to a value of 0.062335 (deg/hr). This is over two orders of magnitude improvement relative to the individual gyro devices, and agrees well with the 0.05 (deg/hr) expected theoretical drift rate.

[0114] In order to form the virtual gyro estimate, it is required to have an estimate of the bias on each of the component gyro devices. When an absolute rate reference is available, the bias estimation can be done individually on a device by device basis. However, such an approach is suboptimal in the sense that it does not take advantage of known correlations between devices. The bias calibration problem will be considered here from an optimal point of view that considers such correlations. As set forth in the appendix an analytic and asymptotic solution to calibration of the Riccati differential equation is possible.

[0115] Many alterations and modifications may be made by those having ordinary skill in the art without departing from the spirit and scope of the invention. Therefore, it must be understood that the illustrated embodiment has been set forth only for the purposes of example and that it should not be taken as limiting the invention as defined by the following claims. For example, notwithstanding the fact that the elements of a claim are set forth below in a certain combination, it must be expressly understood that the invention includes other combinations of fewer, more or different elements, which are disclosed in above even when not initially claimed in such combinations.

[0116] The words used in this specification to describe the invention and its various embodiments are to be understood not only in the sense of their commonly defined meanings, but to include by special definition in this specification structure, material or acts beyond the scope of the commonly defined meanings. Thus if an element can be understood in the context of this specification as including more than one meaning, then its use in a claim must be understood as being generic to all possible meanings supported by the specification and by the word itself.

[0117] The definitions of the words or elements of the following claims are, therefore, defined in this specification to include not only the combination of elements which are literally set forth, but all equivalent structure, material or acts for performing substantially the same function in substantially the same way to obtain substantially the same result. In this sense it is therefore contemplated that an equivalent substitution of two or more elements may be made for any one of the elements in the claims below or that a single element may be substituted for two or more elements in a claim. Although elements may be described above as acting in certain combinations and even initially claimed as such, it is to be expressly understood that one or more elements from a claimed combination can in some cases be excised from the combination and that the claimed combination may be directed to a subcombination or variation of a subcombination.

[0118] Insubstantial changes from the claimed subject matter as viewed by a person with ordinary skill in the art, now known or later devised, are expressly contemplated as being equivalently within the scope of the claims. Therefore, obvious substitutions now or later known to one with ordinary skill in the art are defined to be within the scope of the defined elements.

[0119] The claims are thus to be understood to include what is specifically illustrated and described above, what is conceptionally equivalent, what can be obviously substituted and also what essentially incorporates the essential idea of the invention.

Claims

1. A method of processing signals from multiple rate sensing devices each providing a rate measurement to derive a higher accuracy output signal representing a best estimate of a virtual sensed rate, the method comprising: measuring the rate measurements from multiple rate sensing devices; combining the rate measurements from the multiple rate sensing devices in a recursive processing filter, wherein the weights and smoothing properties of the recursive processing filter are determined to advantage using statistical correlations and noise characteristics of the multiple rate sensing devices.

2. The method of claim 1 wherein combining the rate measurements from the multiple rate sensing devices is performed by a continuous-time steady state Kalman filter to generate a minimum variance rate estimate.

3. The method of claim 1 wherein combining the rate measurements of the multiple rate sensing devices comprises: digitizing the rate measurements of the multiple rate sensing devices; and digitally combining the rate measurements of the multiple rate sensing devices using a digital steady state Kalman filter to generate a minimum variance rate estimate.

4. The method of claim 3 wherein digitally combining the rates of the multiple rates sensing devices using a digital steady state Kalman filter comprises: determining N, wherein N is a number of rate sensing devices being measured; determining Qb and R, wherein Qb and R are N×N covariance matrices whose elements correspond to the correlation factors between the N multiple rate sensing devices; choosing a random walk intensity qω for a true rate, and forming 27Q=[Qb00qω], wherein Q is an (N+1)×(N+1) covariance matrix; generating L=Q1/2HTR−1HQ1/2, wherein L is an (N+1)×(N+1) matrix, wherein H represents an N×(N+1) matrix in a measurement equation y=Hx+v=b+1ω+v; wherein ω represents a true rate along a single axis; wherein b represents a vector of rate biases; wherein x is a state vector x=[bT,ω)]T; wherein v is a vector of angle random walk (ARW) noises that have been stacked into the vector form v=[n11, n12, . . . , n1N]T; and wherein y is a vector representing the measured signals from N gyros; performing an eigenvalue decomposition of L ordered such that its zero eigenvalue is last wherein L=ADAT, and wherein D is a diagonal matrix and A is an orthogonal matrix; extracting a N×N submatrix {overscore (D)} from the matrix D according to 28D=[D_000];extracting a (N×1)×N submatrix {overscore (A)} from the matrix A according to A=[{overscore (A)}v]; generating a steady state Kalman gain as 29K_=S⁡[D_-1/2(I-ⅇ-D_12⁢T)0 00]⁢ST⁢HT⁢R-1=Q1/2⁢A_⁢ ⁢D_-1/2⁡(I-ⅇ-D_1/2⁢T)⁢A_T⁢Q1/2⁢HT⁢R-1 wherein S=Q1/2A and wherein T is a sampling period; and, processing the measurements from the multiple rate sensing devices according to the recursive processing filter, the recursive processing filter having an output and modeled as, 30ξ^⁡((k+1)⁢T)=[ⅇ-D_12⁢T001]⁢ξ^⁡(k⁢ ⁢T)+S-1⁢K_⁢ ⁢y⁡(k⁢ ⁢T)ω^⁡(k⁢ ⁢T)=eN+1T⁢S⁢ ⁢ξ^⁡(k⁢ ⁢T) wherein x represents a vector x=[bT,ω]T; ξ(kT)=S−1x(kT) is a filter state at the kth sampling point; 31eN+1T=[0,…⁢ ,0,1]is a selection vector; and {circumflex over (ω)}(kT) is the output of the recursive filter which represents the best estimate of the virtual sensed rate for the kth sampling point.

5. The method of claim 4 wherein S−1{overscore (K)} is a matrix product having a last row, wherein the last row of the matrix product S−1{overscore (K)} is explicitly set to zero; and wherein

6. The method of claim 2 wherein combining the rate measurements of multiple rate sensing devices performed by the continuous-time steady state Kalman filter comprises: determining N, wherein N is a number of multiple rate sensing devices being measured; determining Qb and R, wherein Qb and R are N×N covariance matrices whose elements correspond to correlation factors between the N multiple rate sensing devices; choosing a random walk intensity qω for a true rate; and forming 33Q=[Qb00qω], wherein Q is a (N+1)×(N+1) covariance matrix; generating L=Q1/2HTR−1HQ1/2 wherein L is a (N+1)×(N+1) matrix, wherein H represents a N×(N+1) matrix in a measurement equation y=Hx+v=b+1ω+v; wherein ω represents a true rate along a single axis; wherein b represents a vector of rate biases; wherein x is a state vector x=[bT,ω]T; wherein v is a vector of ARW noises having been stacked into vector form v=[n11, n12, . . . , n1N]T; and wherein y is a vector representing a measured signals from N gyros; performing an eigenvalue decomposition of L ordered such that its zero eigenvalue is last wherein L=ADAT, wherein D is a diagonal matrix and A is an orthogonal matrix; extracting a N×N submatrix {overscore (D)} from matrix D according to 34D=[D_000]extracting a (N+1)×N submatrix {overscore (A)} from matrix A according to A=[{overscore (A)}v]; generating a steady state Kalman gain K∞ wherein K∞=Q1/2{overscore (A)}{overscore (D)}−1/2{overscore (A)}TQ1/2HTR−1 processing the measurements from the multiple rate sensing devices according to the recursive processing filter having an output, wherein the recursive processing filter is modeled as, 35ξ^.=[-D_12000]⁢ξ^+S-1⁢K∞⁢y ⁢ω^=eN+1T⁢S⁢ ⁢ξ^ wherein S=Q1/2A; wherein x represents the vector x=[bT,ω]T; wherein ξ=S−1x is a filter state; wherein 36eN+1T=[0,…⁢ ,0,1] is a selection vector; and wherein {circumflex over (ω)} is the output of the recursive processing filter which represents the best estimate of the virtual sensed rate.

7. The method of claim 6, wherein S−1K∞ is a matrix product having a last row, wherein the last row of S−1K∞ is explicitly set to zero; and wherein

8. The method of claim 3 wherein combining the rate measurements from the multiple rate sensing devices using a digital steady state Kalman filter to generate a minimum variance rate estimate comprises: sampling N analog signals from N rate sensing devices over a period, T; digitizing the sampled N analog signals into a digital data point for each sampling period T of each rate sensing device; assembling the N digitized signals from the N rate sensing devices into an N-vector; processing the N-vector by multiplying it with a matrix {overscore (K)} operator to generate a {overscore (K)} product signal; processing the {overscore (K)} product signal by multiplying it with a matrix S−1 operator to generate a S−1 product signal; inputting the S−1 product signal during a kth sampling period into a first input of a summing node having two inputs; outputting a difference of the two inputs of the summing node to a delay operator having contents; delaying the contents of the delay operator by one sampling period, T; feeding back a delayed output of the delay operator to a matrix operator modeled by 38&AutoLeftMatch;[ⅇ-D_1/2⁢T001] to generate a feedback signal; inputting the feedback signal to a second one of the inputs to the summing node wherein the feedback signal is added to the S−1 product signal generated in a k+1 sampling period; processing the contents of the delay operator for a kth sampling period by multiplying it with a matrix S operator to generate an S product signal; and extracting a last matrix element of the S product signal as the best estimate of the virtual sensed rate, {circumflex over (ω)}(kT) for the kth sampling point.

9. An apparatus for processing signals from multiple rate sensing devices each providing a rate measurement to derive a higher accuracy output signal representing a best estimate of a virtual sensed rate, comprising: means for measuring the rate from the multiple rate sensing devices; means for combining measured rates from the multiple rate sensing devices in a recursive processing filter, and means for determining the weights and smoothing properties of the recursive processing filter using the statistical correlations and noise characteristics of the multiple rate sensing devices.

10. The apparatus of claim 9 wherein the means for combining the measured rates from the multiple rate sensing devices comprises a continuous-time steady state Kalman filter to generate a minimum variance rate estimate.

11. The apparatus of claim 9 wherein the means for combining the measured rates from the multiple rate sensing devices comprises means for digitizing the rates of multiple rate sensing devices and means for digitally combining the rates from multiple rate sensing devices using a digital steady state Kalman filter.

12. The apparatus of claim 11 wherein the means for digitally combining the rates from the multiple rates sensing devices signals using a digital steady state Kalman filter comprises: means for determining N, wherein N is a number of multiple rate sensing devices being measured; means for determining Qb and R, wherein Qb and R are N×N covariance matrices whose elements correspond to correlation factors between the N multiple rate sensing devices; means for choosing a random walk intensity qω for a true rate, and forming 39Q=[Qh00qω], wherein Q is an (N+1)×(N+1) covariance matrix; means for generating L=Q1/2HTR−1HQ1/2, wherein L is an (N+1)×(N+1) matrix, wherein H represents an N×(N+1) matrix in a measurement equation y=Hx+v=b+1ω+v; wherein ω represents a true rate along a single axis; wherein b represents a vector of rate biases; wherein x is a state vector x=[bT,ω]T; wherein v is a vector of angle random walk (ARW) noises that have been stacked into the vector form v=[n11, n12, . . . , n1N]T; and wherein y is a vector representing the measured signals from N gyros; means for performing an eigenvalue decomposition of L ordered such that its zero eigenvalue is last wherein L=ADAT, and wherein D is a diagonal matrix and A is an orthogonal matrix; means for extracting a N×N submatrix {overscore (D)} from the matrix D according to 40D=[D_000];means for extracting a (N+1)×N submatrix {overscore (A)} from the matrix A according to A=[{overscore (A)}v]; means for generating a steady state Kalman gain as 41K_=S⁡[D_-1/2⁡(I-ⅇ-D_12⁢T)000]⁢ST⁢HT⁢R-1=Q1/2⁢A_⁢ ⁢D_-1/2⁡(I-ⅇ-D12⁢T)⁢A_T⁢Q1/2⁢HT⁢R-1 wherein S=Q1/2A and wherein T is a sampling period; and, means for processing the rates from the multiple rate sensing devices according to the recursive processing filter, the recursive processing filter having an output and modeled as, 42ξ^⁡((k+1)⁢T)=[ⅇ-D_12⁢T001]⁢ξ^⁡(kT)+S-1⁢K_⁢ ⁢y⁡(kT)ω^⁡(kT)=ⅇN+1T⁢S⁢ ⁢ξ^⁡(kT)⁢  wherein x represents a vector x=[bT,ω]T; wherein ξ(kT)=S−1x(kT) is a filter state at the kth sampling point; wherein 43ⅇN+1T=[0,…⁢ ,0,1] is a selection vector; and wherein {circumflex over (ω)}(kT) is the output of the recursive filter which represents the best estimate of the virtual sensed rate for the kth sampling point.

13. The apparatus of claim 12 wherein S−1{overscore (K)} is a matrix product having a last row, wherein the last row of the matrix product S−1{overscore (K)} is explicitly set to zero; and wherein

14. The apparatus of claim 10 wherein the means for combining the measured rates of multiple rate sensing devices comprising a continuous-time steady state Kalman filter further comprises: means for determining N, wherein N is a number of multiple rate sensing devices being measured; means for determining Qb and R, wherein Qb and R are N×N covariance matrices whose elements correspond to correlation factors between the N multiple rate sensing devices; means for choosing a random walk intensity qω for a true rate; and means for forming 45Q=[Qb00qω], wherein Q is a (N+1)×(N+1) covariance matrix; means for generating L=Q1/2HTR−1HQ1/2 wherein L is a (N+1)×(N+1) matrix, wherein H represents a N×(N+1) matrix in a measurement equation y=Hx+v=b+1ω+v; wherein ω represents a true rate along a single axis; wherein b represents a vector of rate biases; wherein x is a state vector x=[bT,ω]T; wherein v is a vector of ARW noises having been stacked into vector form v=[n11, n12, . . . , n1N]T; and wherein y is a vector representing the measured signals from N gyros; means for performing an eigenvalue decomposition of L ordered such that its zero eigenvalue is last wherein L=ADAT, wherein D is a diagonal matrix and A is an orthogonal matrix; means for extracting a N×N submatrix {overscore (D)} from matrix D according to 46D=[D_000]means for extracting a (N+1)×N submatrix {overscore (A)} from matrix A according to A=[{overscore (A)}v]; means for generating a steady state Kalman gain K∞ wherein K∞=Q1/2{overscore (A)}{overscore (D)}−1/2{overscore (A)}TQ1/2HTR−1 means for processing measured rates from the multiple rate sensing devices according to the recursive processing filter having an output, wherein the recursive processing filter is modeled as, 47ξ^.=[-D_12000]⁢ ⁢ξ^+S-1⁢K∞⁢yω^=ⅇN+1T⁢S⁢ ⁢ξ^ wherein S=Q1/2A; wherein x represents a vector x=[bT,ω]T; wherein ξ=S−1x is a filter state; wherein 48ⅇN+1T=[0,…⁢ ,0,1] is a selection vector; and wherein {circumflex over (ω)} is the output of the recursive processing filter which represents the best estimate of the virtual sensed rate.

15. The apparatus of claim 14, wherein S−1K∞ is a matrix product having a last row, wherein the last row of S−1K∞ is explicitly set to zero; and where

16. The apparatus of claim 11 wherein the means for combining the measured rates of the multiple rate sensing devices using a digital steady state Kalman filter to generate a minimum variance rate estimate further comprises: means for sampling N analog signals from N rate sensing devices over a period, T; means for digitizing the sampled N analog signals into a digital data point for each sampling period T of each rate sensing device; means for assembling the N digitized signals from the N rate sensing devices into an N-vector; means for processing the N-vector by multiplying it with a matrix {overscore (K)} operator to generate a {overscore (K)} product signal; means for processing the {overscore (K)} product signal by multiplying it with a matrix S−1 operator to generate a S−1 product signal; means for inputting the S−1 product signal during a kth sampling period into a first input of a summing node having two inputs; means for outputting the difference of the two inputs of the summing node to a delay operator; means for delaying the contents of the delay operator by one sampling period, T; means for feeding back a delayed output of the delay operator to a matrix operator modeled by 50[-D_12001]&AutoLeftMatch; to generate a feedback signal; means for inputting the feedback signal to a second one of the inputs to the summing node such that it is added to the S−1 product signal generated in a k+1 sampling period; means for processing the contents of the delay operator for a kth sampling period by multiplying it with a matrix S operator to generate an S product signal; and means for extracting a last matrix element of the S product signal as the best estimate of the virtual sensed rate, ŵ(kT) for the kth sampling point to generate a minimum variance rate estimate.

17. The apparatus of claim 9 wherein each means comprises a firmware controlled digital processor for performing the recited function.

18. The apparatus of claim 9 wherein each means comprises a programmed portion of a digital software controlled computer for performing the recited function.

19. The apparatus of claim 9 wherein at least one of the means comprises a logic circuit for performing the recited function.

20. The apparatus of claim 9 wherein at least one of the means comprises a firmware controlled digital processor for performing the recited function.

21. The apparatus of claim 9 wherein at least one of the means comprises a programmed portion of a digital software controlled computer for performing the recited function.

22. The method of claim 1 wherein the multiple rate sensing devices are multiple gyros producing signals and the measuring and the combining elements use the signals from the multiple gyros.

23. The apparatus of claim 9 wherein the multiple rate sensing devices are multiple gyros producing signals and wherein the means for measuring and the means for combining use the signals from the multiple gyros.