This is a preprint of a tutorial article.

The final version appears in IEEE Control

Systems magazine. Please cite as:
J. A. Farrell, F. O. Silva, F. Rahman and J. Wendel, "Inertial Measurement Unit Error
Modeling Tutorial: Inertial Navigation System State Estimation with Real-Time Sensor
Calibration," in IEEE Control Systems Magazine, vol. 42, no. 6, pp. 40-66, Dec. 2022,
doi: 10.1109/MCS.2022.3209059.

IMU Error Modeling Tutorial

INS state estimation with real-time sensor calibration

Jay A. Farrell, Felipe O. Silva, Farzana Rahman, J. Wendel

POC: J. Farrell ([email protected])
December 19, 2022

Autonomous vehicle technology is rapidly advancing. Key enabling factors are the
2 advancing capabilities and declining cost of computing and sensing systems that enable sensor
fusion for awareness of the vehicle’s state and surroundings. For control purposes, the vehicle
4 state must be estimated accurately, reliably, at a sufficiently high sample rate, and with a
sufficiently high bandwidth. For systems with high bandwidth, these requirements are often
6 achieved by an aided inertial navigation system (INS) [1], [2], [3], [4], [5], [6]. An INS integrates
data from an inertial measurement unit (IMU) through a kinematic model at the high sampling
8 rate of the IMU to compute the state estimate. An aided INS corrects this state estimate using
data from aiding sensors [for example, vision, LIDAR, RADAR, and Global Navigation Satellite
10 System (GNSS)]. State estimation by sensor fusion may be accomplished by any of a variety of
methods: Kalman filter (KF) [7], [8], [9], [10], [11], extended Kalman filter (EKF) [12], [13],
12 [14], [15], unscented Kalman filter (UKF) [16], [17], [18], particle filter (PF) [19], [20], [21],
and maximum a posteriori (MAP) optimization [22], [23], [24], [25], [26], [27].

14 A data fusion system that combines the IMU and aiding sensor data will be able to
achieve improved performance by real-time calibration if it incorporates an IMU error model in
16 state-space form. The IMU manufacturer supplies a data sheet characterizing the expected IMU
performance. In accordance with the specification standards [28], [29], [30], this performance is
18 typically stated in terms of the Allan Variance (AV). However, it is not immediately clear how
to translate the AV information from such data sheets into a suitable state-space model. Such
20 translation methods have been known and used for several decades [31], [32], [33], [34], [35],
[36], [37], [38], [39], [40]. Despite their importance, a clear tutorial exposition of the underlying
22 ideas, issues, and tradeoffs is not available in the existing literature. Providing such a tutorial
discussion is the purpose of this article. The long history of these ideas and issues relative to
24 successful applications is discussed in the “Aided INS History.”

1
 Application Context

2 The INS design is based on the vehicle kinematic model,

⃗x˙ v (t) = f (⃗xv (t), ⃗u(t)), (1)
where ⃗xv represents the state of the vehicle and ⃗u ∈ ℜ6 represents the system inputs (that is,
4 specific force and angular rate vectors). A typical vehicle state vector might include subvectors
for position, velocity, and attitude. The navigation system numerically solves (1) based on a
6 measurement of the signal ⃗u(t),
˙
⃗xˆv (t) = f (⃗xˆv (t), ⃗uˆ(t)), (2)
where ⃗uˆ(t) is computed from the IMU measurement ⃗u˜(t) using calibration factors that are
8 estimated in real-time. A simplified two-dimensional inertial navigation example illustrating these
ideas is presented in “Simplified INS Example.”

10 For a scalar signal, the model relating the sensor measurement ũ(t) to the desired signal
u(t) is (see Annex B in [28], [29])
ũ(t) = u(t) + d(⃗u(t)) + z(t). (3)
12 The measurement ũ(t) of the desired signal u(t) is corrupted by deterministic errors d(⃗u(t)) and
cumulative stochastic error z(t). Deterministic errors are the sensor imperfections for which an
14 analytical model with unknown deterministic coefficients is sufficient. Some of these coefficients
exhibit only minor variations during the lifetime of the instrument. These can be estimated and
16 compensated in a factory calibration process. Other deterministic errors like turn-on biases can be
estimated in real-time via state augmentation. Forms of deterministic errors might include scale
18 factor error, nonlinearity, g-sensitivity for gyros, nonorthogonal axes, and axes cross-coupling
for vector measurements ⃗u(t) [5], [41], [42]. See the “Deterministic Errors” subsection in the
20 “Discussion of Issues and Tradeoffs.” The focus of this article is the stochastic error denoted by
z(t), which may arise from a variety of physical phenomena (see the “Background” section).
22 The stochastic errors are distinct each time that the instrument is turned on, vary as a function
of time, and cannot be predicted based on the sensor measurement ũ(t).

24 For clarity and simplicity, the majority of this tutorial will treat z(t) as a scalar signal.
The underlying ideas apply to each of the three accelerometers and three gyros in a six-degree-
26 of-freedom IMU for the development of the full IMU stochastic error model.

State Estimation Error Model

28 While the navigation system propagates the vehicle state vector through time by integrating
the nonlinear model of (2), error denoted by δ⃗xv (t) = ⃗xv (t) − ⃗xˆv (t) will develop between the

2
 actual and estimated state vectors. Regardless of the choice of representation of the vehicle
2 attitude (for example, direction cosine matrix or quaternion), the attitude error can be represented
by a vector with three components. Therefore, the vehicle error state contains subvectors for
4 position, velocity, and attitude error, each being a vector with three components such that δ⃗xv (t) ∈
ℜnv with nv = 9. This vehicle state error vector can be estimated in real-time using measurements
6 from the aiding sensors [10], [13], [12], [31], [43], [44]. The estimation algorithm incorporates
a linearized state-space model for the error state,

δ⃗x˙ v (t) = F (t)δ⃗xv (t) + G(t)δ⃗u(t), (4)

8 where F (t) = ∂f (⃗
∂⃗
x,⃗
x
u)
, G(t) = ∂f (⃗
∂⃗
x,⃗
u
u)
, and δ⃗u(t) = ⃗u(t) − ⃗uˆ(t). See the example
ˆ(t), ⃗
⃗
x ˆ(t)
u ˆ(t), ⃗
⃗
x ˆ(t)
u
F (t) and G(t) matrices in “Simplified INS Example.” This state-space error model is not
10 complete until the IMU error model for δ⃗u(t) is specified. For state estimation, it is required
that the IMU error model be in state-space form.

12 The real-time state estimation process is designed to estimate the augmented state vector
h i⊤
⊤ ⊤
⃗x(t) = δ⃗xv (t) ⃗xd (t) ⃗xz (t) ⊤ ∈ ℜnx , (5)

comprised of the vehicle error state vector δ⃗xv (t) ∈ ℜnv , the vector ⃗xd (t) ∈ ℜnd augmented to
14 enable calibration of the IMU deterministic errors, and the vector ⃗xz (t) ∈ ℜnz augmented to
enable calibration of the IMU stochastic errors. The total number of error states nx = nv +nd +nz ,
16 where nv is the number of vehicle error states, nd is the number of states augmented to calibrate
deterministic errors, and nz is the number of states augmented to calibrate stochastic errors.
18 The process of state augmentation as it relates to the problems of interest is discussed in “State
Augmentation.” The presentation herein will focus entirely on the stochastic IMU errors whose
20 cumulative effect in (3) is denoted by ⃗z(t). The majority of this article will focus on a single IMU
output for which a scalar z(t) is sufficient. The “Discussion of Issues and Tradeoffs” section
22 considers the vector ⃗z(t) case.

Denote the state-space model for scalar z(t),

⃗x˙ z (t) = Az ⃗xz (t) + Bz ω

⃗ z (t), (6)
z(t) = Cz ⃗xz (t) + ηz (t), (7)

24 where Az ∈ ℜnz ×nz , Bz ∈ ℜnz ×p , and Cz ∈ ℜ1×nz . The parameter p represents the number of
distinct and independent noise processes in the differential equation portion of the the IMU error
26 model. The parameter nz represents the number of states in the IMU stochastic error model. The
random signals ω ⃗ z (t) and ηz (t) are mutually independent Gaussian white noise processes with
28 power spectral densities (PSD) Sωz ∈ ℜp×p and Sηz ∈ ℜ, respectively. The elements of ω ⃗ z (t)
are assumed to be independent, which yields Sωz being diagonal.

3
 The designer must be judicious in the choice of the model structure (particularly nz and
2 p) as the overall model will have 6 nz states and 6 p independent noise sources.

Power Spectral Density for Linear State-Space Systems

4 Corresponding to the state-space model in (6)-(7), the frequency domain model is

Z(s) = T (s) Ωz (s) + ηz (s), (8)

where s is the Laplace variable. The transfer function model from ωz (t) to z(t) is

T (s) = Cz (sI − Az )−1 Bz ,

which has one row and p columns. The symbols Z(s), Ωz (s), and ηz (s) represent the Laplace
6 transforms of the signals z(t), ωz (t), and ηz (t). Therefore, the PSD corresponding to signal z(t)
is
Sz (ω) = T (jω) Sωz T (−jω)T + Sηz . (9)

8 Assuming that all elements of the driving noise vector ωz (t) and the output noise ηz (t) are
mutually independent and white, this simplifies to
p
X
Sz (ω) = Ti (jω) Ti (−jω) Sωzi + Sηz , (10)
i=1

10 with
Ti (s) = Cz (sI − Az )−1 Bzi , (11)

where Ti (s) is the (scalar) transfer function from the i-th component of ω ⃗ z (t) to z(t) and
nz ×1
12 Bzi ∈ ℜ is the i-th column of Bz . Each Ti (s) is the ratio of a numerator and denominator
polynomial in s, and Sz is a positive real function of ω. Therefore, each Ti (jω)Ti (−jω) and
14 Sz (ω) will always be the ratio of polynomial functions of only the even powers of the Laplace
variable s. Examples demonstrating this fact are presented in “Finite Dimensional Linear State-
16 Space Systems have Even Power Spectra.”

The fact that Sz is a positive real function leads to one of the main challenges in the IMU
18 error modeling approach, because some of the IMU stochastic error components have PSDs that
cannot be exactly fit by the terms in the summation of (10). Therefore, the designer must make
20 judicious choices in the approximate state-space model to achieve satisfactory tradeoffs. This
will be further discussed in the section entitled “Modeling via Independent Noise Sources.”

4
 Available Error Specification Information

2 To characterize the quality of an IMU (per the IEEE specifications [28], [29], [30]) the
manufacturer provides the AV plot, Allan Standard Deviation (ASD) plot, or parameters extracted
4 from them. An ASD graph can both help instrument designers understand and improve their
sensors and communicate expected performance to perspective users. Within the context of this
6 article, the main topic is how an INS designer can use information from an ASD plot to specify
the parameters of the state-space model in (6)-(7).

8 Example ASD plots are shown in Figures 1 and 2. Figure 1 displays the ASD plots for
three gyros in a Crossbow µNav IMU [45]. These ASD are computed from data provided by the
10 authors of [36]. The blue, black, and green asterisks (‘*’) mark the ASD data points. The ASD
plot in Figure 2 is computed from manufacturer-supplied data obtained from an IMU mounted
12 on a large marble slab on top of a vibration isolation system. Each blue ‘x’ in Figure 2 marks
an ASD data point. The horizontal axis of each ASD plot is the cluster time (or size), with
14 the symbol τ measured in seconds. Note that the ASD plots in these two figures have both
similarities and differences. Both decrease with the characteristic slope of -1/2 for small cluster
16 sizes, as indicated by the red dashed tangent line in each figure. Then, both level off to a slope
of zero at values indicated by the dashed cyan tangent lines. For larger cluster sizes, the ASD in
18 Figure 2 increases with a slope of 1/2, as indicated by the tangent line drawn with black dashes.
For cluster times as large as τ = 1000, the ASD plots in Figure 1 do not (strongly) exhibit this
20 increase with slope 1/2. The value of τ and the ASD value at which these changes occur is
distinct for each instrument. They specify certain parameters that are useful both for comparing
22 the performance of inertial instruments and for evaluating tradeoffs in the construction of the
IMU stochastic error model.

24 Summary

This tutorial discusses issues and tradeoffs related to an example navigation system design
26 method: (1) using the AV information to specify the continuous-time parameters (for example, p,
nz , Az , Bz , Cz , Sz , and Sη ) for an IMU state-space stochastic error model; (2) transforming this
28 continuous-time model to the discrete-time, state-space error model parameters (Φ, Qzd , H, Qηd )
that are required for implementation of the state estimator; and, (3) verifying the IMU model
30 relative to the AV information. Some previous articles have addressed some of the aforementioned
topics [32], [34], [36], [39], [43], [44], [46], [47]. The goal of this article is to clearly and
32 comprehensively present the background and main ideas in a tutorial fashion, using notation and
terminology consistent with instrument specification standards [28], [29]. Examples are included
34 throughout to clarify issues.

5
 Problem Statement

2 The purpose of this article is to discuss methods, issues, and tradeoffs related to developing
a model that quantitatively communicates to the state estimation algorithm the nature of the
4 stochastic portion of the IMU error. The inputs to this error model will be random signals.

The available information from the manufacturer for the model development is the AV
6 or ASD characterization. The actual IMU output data from the experiments that produce these
characterizations are not available to the designer; therefore, system identification methods are
8 not applicable.

Because the error state estimation algorithms are formulated in state-space form, the IMU
10 error models will also be in state-space form. The inputs are modeled as independent Gaussian
white noise processes. The challenge is to construct these models such that they have the same
12 output statistical characteristics (that is, AV) as the IMU.

Note that such stochastic state-space models are not unique. In fact, this article is not
14 intended to propose a particular model; although a specific ASD plot is discussed and a model
is given as a tutorial example. Instead, the goal of this article is to present clearly the approach
16 that is used or hinted at (with various specific models) in various articles, books [9], and standards
[30], [28], [29]. The authors view the many flavors of this approach as being industry-standard.
18 Unfortunately, most publications describing the method are not publicly available.

From a high-level, an outline of the method is as follows: (1) A continuous-time state-

20 space model is constructed using information from the AV/ASD plot. (2) The continuous-time
IMU stochastic model is transformed to an equivalent discrete-time model. (3) The discrete-
22 time model is used in simulation to produce data, from which an ASD plot is computed for
comparison with the instrument’s ASD plot. (4) When the designer is satisfied with the IMU
24 error model, it is appended to the vehicle state error model and used for the design of the INS
error state estimator.

26 Background

Because the manufacturer supplied AV/ASD information is the starting point for the
28 stochastic error model development, this section reviews the AV and its relationship to the
power spectral density. A brief history of the AV is included in “A Brief Historical Review of
30 the Allan Variance.”

6
 Allan Variance

2 The AV is a well-known time domain analysis technique that was originally developed to
characterize and study the frequency stability of oscillators [48], [49], [50]. Due to its relative
4 simplicity, it has been successfully adopted to communicate IMU performance specification and
to characterize their stochastic errors [28], [29], [32], [34], [36], [39], [43], [44], [46], [47].

6 Given a set of data, the process for computing the AV is as follows. Let D = {ũi }Li=1
be a (detrended) set of specific force (or angular rate) data, measured at a constant sampling
8 interval T for a stationary (that is, motion isolated) IMU. For each n ∈ [1, L/2], the AV is
computed for values of the cluster time τ = nT ranging from T to L T /2. For a given n, at
10 each time instant ti ∈ [T, 2T, . . . , (L − n) T ], a group of n consecutive data points (beginning
at ti ) form a cluster: {ũj }i+n−1
j=i . The average value is computed for each such n-point cluster,
1
Pn−1
12 ūi (τ ) = n j=0 ũi+j . The AV for duration τ is then computed as the average of the (L − 2n)
squared cluster differences [28], [49]:
L−2n
1 X
σ̂u2 (τ ) = [ūi+n (τ ) − ūi (τ )]2 . (12)
2(L − 2n) i=1
14

Since some IMUs (especially those that are high-grade) provide the integral of specific
16 force (or angular rate), denoted as θ̃i , ūi (τ ) may alternatively be defined as

θ̃i+n − θ̃i
ūi (τ ) = . (13)
τ
Substitution of (13) into (12) yields
L−2n
1 X
σ̂u2 (τ ) = 2 (θ̃i+2n − 2θ̃i+n + θ̃i )2 , (14)
2τ (L − 2n) i=1

18 which is an alternative formula for computing the AV [28].

For graphical analysis, the square root of the AV, σ̂u (τ ), called the ASD, is typically plotted
20 on a log-log scale with cluster time τ along the horizontal axis. Due to the finite length of dataset
D, the number of clusters with duration τ will decrease as τ increases; therefore, the standard
22 deviation of the computed ASD, σ̂u (τ ), increases with n (or τ ) as [51], [52]
r
n
σ [σ̂u (τ )] = κ σ̂u (τ ), (15)
L
√
where κ is an empirical constant that is generally approximated as κ ≈ 1/ 2 for IMU error
24 analysis [53], [28], [29], [34], [54].

7
 Power Spectral Density and Allan Variance

2 The AV is related to the two-sided PSD by

Z ∞
2 sin4 (π f τ )
σu (τ ) = 4 Su (f ) df. (16)
0 (π f τ )2
The text following (C.1) in IEEE Standard 952-1997 [28] interprets this equation as the AV being
4 proportional to the total noise power in the signal u when passed through a transfer function that
is determined by the method that is used to create and operate on the clusters. The derivation
6 of (16) can be found on p. 79 in [55]. There is no inversion formula for (16) (see [53]). In this
√
expression, Su (f ) = Su (s)|s=j2πf , where s ∈ C is the Laplace variable, j = −1, and f ∈ ℜ
8 has units of Hertz.

Modeling via Independent Noise Sources

10 When the power spectrum is represented as a power series in frequency f , it has the form

B2
2 K2
Su (f ) = · · · + N + + + ··· . (17)
2πf (2πf )2
This form of the PSD is convenient. By the principle of superposition, it corresponds to the
power spectrum of the signal

u(t) = · · · + zN (t) + zB (t) + zK (t) + · · · (18)

where the signals zN (t), zB (t), and zK (t), are mutually independent, zero mean, noise processes.
With this assumption, applying (16) to (17) yields an AV having the form

σu2 (τ ) = · · · + σz2N (τ ) + σz2B (τ ) + σz2K (τ ) + · · · , (19)

12 where the specific functional form of each AV term can be computed and is available in various
sources [28], [29], [30], [34], [55]. The functional form of each AV term is easily associated
14 with a portion of the ASD graph.

The “Continuous-Time State-Space Models” section describes the method for making this
16 association and defining a continuous-time state-space model for each term. The state-space
model for each term is driven by its own independent, Gaussian, white, driving noise resulting
18 in each of the signals zN (t), zB (t), and zK (t) being mutually independent. These state-space
models can be exact for the terms that correspond to even functions of f in (17). However
20 (as previously stated and as exemplified in a sidebar) the power spectrum terms that are odd
B2
functions of the frequency f (for example the term 2πf ) cannot be exactly modeled by any finite
22 dimensional, linear, state-space model. Therefore, these terms must be approximately modeled,
carefully balancing tradeoffs that are discussed later.

8
 Any number of terms may be included in the power series representation of (17). This
2 results in the same number of terms in the signal model of (18) and the AV model of (19).
Each term represents a different type of noise coming from an independent source. The typical
4 shape of the ASD graph is depicted in Figure 3 with five independent noise sources (see also
Figure C.8 in [28]). In the ASD plot, each noise type is associated with a characteristic slope
6 that facilitates identification of that noise type and its model parameters. Not all noise types are
evident in each instrument. When present, the model parameters and range of τ over which the
8 noise term is dominant may be different for each instrument.

The N , B, and K terms are typically dominant in commercial-grade IMU’s (see for
10 example Figures 1 and 2). Instrument design choices (for example quantization approach and
sample period) cause the stochastic error to appear as white noise for small τ . This white noise
12 is accounted for in the random walk noise term (that is N ). However, the stochastic error is not
truly white. As the cluster time τ increases, the ASD plot may exhibit bias instability (B), rate
14 random walk (K), and other noise types. For the ASD plot to exhibit these other noise types,
the IMU data set used to generate the ASD plot must be very long. When the INS is designed to
16 work with aiding measurements that are expected to occur frequently (for example, several times
per minute), the state estimator will have the aiding information that it needs to maintain the
18 INS calibration in real-time, while those aiding measurements are available. The ASD plot out
to several minutes (for example, hundreds of seconds) is of interest for analyzing performance
20 during intervals when the aiding measurements are not available. However, the specific shape of
the ASD curve for very large τ is typically uncertain and not of interest.

22 Columns 1 and 2 of Table 1 include the specific names of the N , B, and K noise terms for
gyros and accelerometers. Columns 3 and 4 of Table 1 summarize the relationships between the
24 AV and PSD for these noise types, as derived in [28], [34]. The N , B, and K terms will be the
focus of the discussion in the “Continuous-Time State-Space Models” section. The underlying
26 ideas of the approach extend to other types of noise when the ASD for a particular instrument
exhibits them.

28 Continuous-Time State-Space Models

This section considers the development of continuous-time state-space models that approx-
30 imately reproduce the ASD plots and the PSD of (17). The ideas throughout will be illustrated
using the example ASD in Figure 2. Figure 1 will only be discussed in reference to the B and
32 K power series terms. The overall model will have the form,

z(t) = zN (t) + zB (t) + zK (t), (20)

9
 where zN (t), zB (t), and zK (t) are the IMU stochastic error signals associated with coefficients
2 N , B, and K, respectively.

Random Walk Errors: Angular and Velocity: zN (t)

4 The PSD term N 2 in (17) is constant with respect to frequency f , which corresponds to
the power spectrum of white noise [56]. Therefore,

zN (t) = ωN (t), (21)

6 where ωN (t) is white Gaussian random noise with PSD

SN = N 2 . (22)

In the literature and manufacturer specifications, this type of error is called angular random walk
8 error for gyros and velocity random walk error for accelerometers.

Applying the transformation in (16) to SzN (f ) = N 2 yields [28],

N2 N
σz2N (τ ) = or σzN (τ ) = 1/2 , (23)
τ τ
10 which is summarized in the corresponding row of Table 1. This shows that, on an ASD plot,
angular/velocity random walk will be represented by a line with a slope of − 12 , as shown in
12 Figure 4.

The value of the random walk parameter N can be approximately determined from the
14 manufacturer supplied ASD plot. This is accomplished by identifying the range of τ on the ASD
plot that has a slope of − 12 and drawing a line tangent to it. In Figure 2, a red dashed tangent
16 line is drawn for τ ∈ [0.01, 30]. From (23), it is clear that σzN (τ ) τ =1 = N . Therefore, the value
of N can be extracted from the ASD plot as the value of the tangent line (with slope of − 21 ) at
18 τ = 1s. For the example of Figure 2, the result is N ≈ 0.0033 m/s3/2 .

Random Walk Errors: Rate and Acceleration: zK (t)

2
K
20 The term (2πf )2
= Ks sK∗ s=j2πf in (17) corresponds to a linear system with transfer function
T (s) = 1s . A state-space model is
żK (t) = ωK (t), (24)

22 where zK (t) is the output and the input ωK (t) is white Gaussian noise with PSD

SK = K 2 . (25)

In the literature and manufacturer specifications, this type of error is called rate random walk
24 error for gyros and acceleration random walk error for accelerometers.

10
 Given (24) and (25), the PSD of zK (t) is
K2
 
SzK (f ) = T (s)T (s∗ ) K2 = , (26)
s=j2πf (2πf )2
2 which has the desired form corresponding to the third term in (17). Using (16) on this SzK (f )
yields [28] r
2
K τ τ
σz2K (τ ) = or σzK (τ ) = K , (27)
3 3
4 which is summarized in the corresponding row of Table 1. Equation (27) shows that on an ASD
plot, the rate/acceleration random walk error will be represented by a line with a slope of + 12 ,
6 as shown in Figure 5.

The rate/acceleration random walk parameter K can be approximately determined from

8 the manufacturer-supplied ASD plot. The first step is to identify the range of τ on the ASD plot
that has a slope of + 21 (if it exists) and draw a line tangent to it. In Figure 2, the dashed black
10 line drawn for τ ∈ [3, 500] seconds, has slope + 12 and is approximately tangent to the ASD curve
at values of τ ≥ 100 seconds. Because τ is large, this portion of the ASD plot usually has a
12 higher degree of uncertainty [as discussed relative to (15)]. The second step uses the tangent line
to estimate K. From (27), it is clear that σzK (τ ) τ =3 = K. Therefore, an easy way to estimate
14 the value of K from the ASD plot is to find the value of the straight line approximation when
τ = 3 s. In the example of Figure 2, depending on how the analyst determines the straight-line
16 approximation, K ≈ 0.00014 m/s5/2 .

Based on the ASD plots in Figure 1, the gyros in the µNav unit may not require inclusion
18 of angular rate random walk noise for cluster times up to 1000 seconds.

Cumulative Error Model: N, K

20 Since the angular random walk and rate random walk errors (or velocity random walk
and acceleration random walk errors) each have even power spectra, it was straightforward to
22 establish state-space models to reproduce the corresponding terms in the power spectrum and
their portions of the ASD plot.

24 Based on the two previous sections, the state-space model would be

żK (t) = ωK (t) (28)

zN K (t) = zN (t) + zK (t), (29)

where zN (t) = ωN (t) is white random noise with PSD N 2 and ωK (t) is a white random noise
26 with PSD K 2 . The random signals ωN (t) and ωK (t) are independent, which results in zN (t) and

11
 zK (t) being independent. The ASD for this model is

2 N 2 K 2τ
σN K (τ ) = + . (30)
τ 3
2 The ASD plot for this model using the values N = 0.0033 m/s3/2 and K =0.00014 m/s5/2 is
shown in Figure 2 as the solid green line.

4 Bias Instability: zB (t)

Some ASD plots (such as those in Figure 1) do not exhibit the +1/2 slope associated with
6 zK (t) but do have a wide flat region for larger values of τ . This flat region cannot be well
modeled by either the N or K terms. For other instruments, the ASD plot may exhibit a wide
8 and flat portion between the values of τ corresponding to the N and K regions. In this case, the
AV (and ASD) in (30) for the N K-model may be too small in this middle range of τ . In either
10 of these circumstances, there is enough bias instability such that performance may be improved
by accounting for it in the model.
B2
The error term zB (t) corresponding to SzB (f ) = 2πf is generally referred to as the bias
B2
instability (or flicker noise) [28], [29], [34], [44], [39]. Applying (16) to SzB (f ) = 2πf for f ≤ f0
(and 0 otherwise) yields [28]

2 B2 h sin3 (x) i
σz2B (τ ) = ln(2) − (sin x + 4x cos x) + Ci(2x) − Ci(4x) , (31)
π 2x2
12 where x = πf0 τ, Ci is the cosine integral function [57], and the parameter f0 is defined as the
cut-off frequency [28].

14 The bias instability ASD plot is shown in Figure 6. The figure shows that σzB (τ ) grows
for small τ reaching a plateau for τ > f10 . Therefore, the value of τ ≈ f10 defines the portion of
16 the ASD plot for which the bias instability (or flicker noise) contributes its maximum value to
the SD plot. In this region, it is possible to show that the sine and cosine terms in (31) approach
18 zero, so that in the flat region,
2 B 2 ln(2)
σz2B (τ ) ≈ or σzB (τ ) ≈ 0.664 B. (32)
π
These equations provide a simple method for extracting an approximate value of B from the
20 ASD plot. In this approach (see for example, p. 6 in [35], p. 21 in [44], p. 10 in [54], and p.
114 in [55]), as can be inferred from Section B.4.5 in [28], the value of B would be selected so
2
22 that 2 B πln(2) approximates the plot of σz2 (τ ) for the values of τ for which the ASD plot is flat.
For the ASD plot in Figure 1, the cyan horizontal lines approximate the minimum ASD values
24 of 9.5e-3 and 1.40e-2 deg/s, which correspond to values of B between 1.43e-2 and 2.11e-2

12
 deg/s. For the ASD plot in Figure 2, the minimum ASD value of 7.4e−4 m/s2 corresponds to
2 B = 1.11e−3 m/s2 .
2B
Because the power spectrum of the bias instability term (that is 2πf ) is not an even power of
4 s = j2πf , there is no finite-order linear state-space model that fits it exactly. As a consequence,
the navigation system designer must select a state-space model to approximate the bias instability
6 error effects. This is somewhat of an art, as each IMU is distinct, each application has different
specifications, and each designer may have different ideas about suitable models and tradeoffs.

8 Various methods have been suggested to approximately account for the bias instability.
These include first-order Gauss-Markov [4], [32], [36], [39], [44], [58] and higher-order
10 autoregressive models [43], [46], [47]. One important tradeoff is that as the dimension of
the state-space model increases, the fidelity of the approximation may increase, but so does
12 the required real-time computational load of the state estimation algorithm. In addition, more
elaborate models may not be robust to unmodeled dynamics and nonlinearity, especially when
14 some added states are weakly observable. These topics are analyzed further in the “Discussion
of Issues and Tradeoffs” section.

16 To exemplify the idea, the next section considers a first-order Gauss-Markov model, which
uses exponentially correlated noise to model the bias instability error.

18 Gauss-Markov Error Model

A first-order continuous-time Gauss-Markov model is [9], [56]

żG (t) = −µB zG (t) + ωB (t), (33)
20 with
1
µB = , where TB > 0. (34)
TB
The symbol TB represents the correlation time of the process. The symbol ωB (t) represents a
22 white driving noise with PSD SB .

The transfer function corresponding to (33) is

1
T (s) = ,
s + µB
which yields the PSD
SB
SzG (ω) = 2 .
ω + µ2B
Applying (16) to SzG (s) yields [28]
SB TB2
 
2 TB  − Tτ − T2 τ
σzG (τ ) = 1− 3 − 4e B + e B . (35)
τ 2τ

13
 A plot of σzG (τ ) is shown in Figure 7. Some special cases are noteworthy.

2 • For smaller cluster times, where τ << TB ,

SB τ
σz2G (τ ) ≈ , (36)
3
so that the ASD plot has a slope of + 21 for small τ .
4 • When τ = 1.89 TB , the curve is flat with
p
2
σz2G (1.8 9TB ) = 0.4365 SB TB (37)
√
or σzG (1.89 TB ) = 0.4365 SB TB .
6 • For larger cluster times, where τ >> TB ,
SB TB2
σz2G (τ ) ≈ , (38)
τ
so that the ASD plot has a slope of − 21 for large τ .

8 The first-order scalar Gauss-Markov process can be used (approximately) to model the flat portion
(that is, bias instability) of the ASD plot.

10 If the manufacturer only provides the values of B and TB , then the value of µB can be
computed using (34). The value of SB is selected by setting σz2B from (32) equal to σz2G [as
12 given by (37)], and solving for SB :
2 B 2 ln(2)
SB = . (39)
π(0.4365)2 TB
With µB and SB known, the state-space model of (33) is completely specified.

14 If instead, the manufacturer provides the ASD plot, and the bias instability is significant
enough to warrant inclusion in the model, then the analyst can first select TB so that 1.89 TB
16 lies near the flat portion of the ASD plot. Then, choose SB so that the value of σzG (1.89 TB ), as
defined in (37), approximates the value of the ASD plot in its flat region. For the ASD plot in
18 Figure 2, the minimum value is 7.4e−4 m/s2 at τ = 60 s. These values correspond to TB = 31.7
s, SB = 9.0e−8 m2 /s5 , and B = 1.11e−3 m/s2 .

20 Cumulative Error Model: N , B, K

Consider the two-state state-space model structure, where z(t) is modeled by (20), zN (t)
22 is modeled by (21), zK (t) is modeled by (24), and the bias instability term in (20) is modeled
by zB (t) = zG (t) as defined in (33). This yields a two-state model in the form of (6)-(7) with
" # " #
−µB 0 1 0 h i
Az = , Bz = , Cz = 1 1 , (40)
0 0 0 1

14
 h i⊤ h i⊤
where ⃗xz (t) = zG (t), zK (t) , ω ⃗ z (t) = ωB (t), ωK (t) , and ηz (t) = ωN (t). The
2 continuous-time process noise PSD matrix Sωz of ω⃗z is
" #
SB 0
Sωz = (41)
0 SK

and the measurement noise PSD is Sηz = SN .

4 The ASD for this state-space model is

1/2
σz (τ ) = σz2N (τ ) + σz2G (τ ) + σz2K (τ ) , (42)

where the three ASD terms on the right are defined in (23), (27), and (35).

6 Figure 8 builds on the ASD plot in Figure 2. The data (blue x’s) and K and N tangent
lines are the same. Figure 8 also shows the ASD of (42) using the two sets of parameters, as
8 summarized in the second and third rows of Table 2 (that is, Untuned and Manually Tuned).

The untuned parameters from the second row of the table are those stated at the end of
10 the “Gauss-Markov Error-Model” section. These values result in the green dashed curve, which
is clearly too high in the region of the curve near its minimum. The suggested approaches stated
12 above for choosing N , B, and K independently neglect the fact that for each τ , the AV for each
type of noise is additive:

σz2 (τ ) = σz2N (τ ) + σz2B (τ ) + σz2K (τ ). (43)

14 Therefore, the parameters N , B, K, and TB should instead be adjusted jointly so that the
three term AV (or ASD) model fits the ASD plot for the instrument data. Manually adjusting
16 the values of TB and B to those in the third row of Table 2 results in the solid green curve,
which is a better fit. An optimization-based approach to selecting the parameters is discussed in
18 “Optimization-based Parameter Selection for a Given Model.”

Figure 8 includes the curve for the Gauss-Markov ASD of (35) plotted as a cyan dashed
20 line, using the tuned values of TB and SB .

Summary

22 Many alternative choices for the structure of the state-space model exist: no states, just
zN ; a single state zG and zN ; a single state zK and zN ; generalizations of the above two-state
24 model; and, higher dimensional models. Various topics related to selection of the structure of
the state-space model are presented in the “Discussion of Issues and Tradeoffs” section.

15
 Discrete-time Equivalent Model

2 The previous section discussed the process of determining a continuous-time state-space

model with the form of eqns. (6)-(7) to approximate the IMU error characteristics, as quantified
4 by the ASD (or power spectrum) plot. Since the INS and EKF are implemented in discrete-time,
the state-space IMU error model must be transformed to the equivalent discrete-time form:

⃗xz (k + 1) = Φ ⃗xz (k) + ω

⃗ z (k), (44)
z(k) = H ⃗xz (k) + η(k), (45)

6 where discrete-time k index corresponds to continuous-time t = k T , ω ⃗ z (k) ∼ N (0, Qzd ) is a

white Gaussian random variable with covariance Qzd , and η(k) ∼ N (0, Qηd ) is a white Gaussian
8 random variable with covariance Qηd . The processes ω ⃗ z (k) and η(k) are independent. Within
this article, equivalent means that the continuous-time and discrete-time stochastic error models
10 produce the same first- and second-order statistics at the IMU sampling times.

The following subsections describe how to compute the discrete-time model parameters
12 (Φ, Qzd , H, Qηd ) that are required for the discrete-time estimator implementation from the
continuous-time model parameters (Az , Bz , Sωz , Sηz , Cz ).

14 Discrete-Time Equivalent to Equation (6)

The process of transforming the continuous-time model of (6) to the discrete-time state-
16 space model of (44) is described in various references (see for example Section III.D in [59] or
Section 4.7 of [3]) with

Φ = eAz T and (46)

ZT
Qzd = Φ(T, s) Bz Sωz (s) BzT Φ(T, s)T ds, (47)
0


18 where T = tk − tk−1 is the IMU sampling interval and Φ(T, s) = exp Az (T − s) . Both Φ
and Qzd can be computed simultaneously using a method by Van Loan [60], as described in
20 Appendix I of [59].

Depending on the method used to model the bias instability and the state definition, the Az
22 matrix in the continuous-time state-space model may have various forms. For the Cumulative
h i ⊤
N -B-K Error Model with the discrete-time state vector, ⃗xz (k) = zG (k), zK (k) , the Az
24 matrix in (40) transforms to " #
e−µB T 0
ϕ= . (48)
0 1

16
 Since Az and Sωz are constant and diagonal, and Bz is the identity, (47) simplifies to
" RT
#
SB 0 exp − 2µB (T − p) dp 0
Qzd =
0 SK T
" #
SB
1 − e−2µB T


2 µB
0
= . (49)
0 SK T
1
2 Because T ≪ TB and µB = TB
,it is clear after expanding the exponential term that
" #
SB T 0
Qzd ≈ Sωz T = . (50)
0 SK T

Discrete-Time Equivalent to Equation (7)

4 The objective of this section is to determine an equivalent discrete-time measurement model

having the form of (44)-(45) that is equivalent to (6)-(7).

6 Because the method of the previous section was designed such that (6) and (44) produce
the same first- and second-order statistics for the state, it is the case that H = Cz . The only item
8 left for consideration is to determine how to compute the covariance Qηd from the PSD SN so
that the overall effect (in terms of the first two moments) of z(k) on the integrated state in the
10 augmented state-space model is the same as that of z(t) at the sampling periods.

Using the principle of superposition, assume that the driving noise terms ωB (t) and ωK (t)
12 [and therefore ωz (t)] are zero, and the covariance of the initial conditions is zero; so that the
only remaining nonzero error is ηz (t), which equals ωN (t) = zN (t) [according to the model in
14 (21) and (40)]. Consider the simple case where the kinematic model of (1) is

ẋv (t) = u(t). (51)

This could, for example, correspond to u being an angular rate (or acceleration) measurement
16 and xv representing the angle (or velocity). In continuous-time, the navigation system would use
the measurement ũ(t) to compute x̂v (t),

x̂˙ v (t) = ũ(t). (52)

18 With the assumptions stated at the beginning of the paragraph,

ũ(t) = u(t) + ηz (t). (53)

Define the error signal, e(t) = xv (t) − x̂v (t). Based on (51)-(53),

ė(t) = −ηz (t), (54)

17
 where according to (21)-(22) the PSD of ηz (t) = ωN (t) is SN . In this special case, e(t) is a
2 continuous-time random walk process. Due to the assumption that the initial covariance of e(t)
is zero [that is, cov(e(0)) = Pe (0) = 0], the covariance function [that is, cov(e(t)) = Pe (t)] is

Pe (t) = SN t, for any t ≥ 0. (55)

4 This result is derived in many textbook discussions of the continuous-time random walk
processes, for example the discussion of (4.85) in [3].

6 The discrete-time model that is equivalent to (51) is

xv (k + 1) = xv (k) + u(k) T. (56)

The discrete-time model that is equivalent to (52) is

x̂v (k + 1) = x̂v (k) + ũ(k) T. (57)

8 Given the discrete-time measurement model,

ũ(k) = u(k) + η(k), (58)

where cov(η(k)) = Qηd [as defined in (45)], the error signal e(k) = xv (k) − x̂v (k) has the time
10 propagation model
e(k + 1) = e(k) − η(k) T. (59)

In this special case, e(k) is a discrete-time random walk process. Equation (59) allows
12 computation of the discrete-time error covariance caused by η(k) [that is, Pe (k) = cov(e(k))]
as

Pe (k + 1) = Pe (k) + T 2 Qηd , for any k ≥ 0, (60)

14 where Qηd is defined in (45). Due to the assumption that the initial covariance of e(k) is zero,

Pe (k) = k T 2 Qηd , for any k ≥ 0. (61)

Given the equivalence objective stated in the first sentence of this section, the continuous
16 and discrete time must result in the same values for the error covariance at the discrete sample
times. Equating the covariance of the continuous-time and discrete-time random walk error
18 processes at the sampling times – that is, setting (55) equal to (61) – yields,

Pe (t) t=kT
= Pe (k)
SN k T = k T 2 Qηd ,

which provides the equation (see also “Discussion of Equation (62)” sidebar):
SN
Qηd = . (62)
T

18
 This equation relates the covariance Qηd of the discrete-time IMU measurement noise (needed
2 for the state estimator design) to the PSD SN of the continuous-time measurement noise (derived
from the ASD plot). Substituting (62) into (61) shows that the covariance of the integrated error
4 signal (e(k)) increases linearly with time in proportion to the PSD SN , as is expected for a
random walk.

6 Together, (46)-(47), (62), and H = Cz are the conversions needed to transform the
continuous-time state-space error model to its equivalent discrete-time error model as necessary
8 for the EKF design.

ASD Verification of the State-Space Model

10 The previous sections have presented methods to develop continuous and discrete-time
state-space models to approximate the IMU stochastic errors, as characterized by the AV method.
12 If the method that a designer uses is valid, then data generated by the discrete-time state-space
model should result in an ASD plot that closely approximates the ASD plot provided by the
14 IMU manufacturer. This section contains an example to demonstrate this verification process.

The example starts from the ASD plot of Figure 2. The model, parameter selection, and
16 data generation methods are as follows:

1) The AV parameters are extracted from the ASD plot using the method described in the
18 “Continuous-Time State-Space Models” section. The parameters for this example are stated
in the third row of Table 2 (that is, manually tuned).
2) The value of µB = T1B = 0.05 s−1 by (34), SB = 1.8528e−8 m2 /s5 by (39), and SK =
1.9600e−8 m2 /s5 by (25). By (41),
" #
1.8528e−8 0.00
Sωz = .
0.00 1.9600e−8

For Az , Bz , and Cz defined in (40), using (48) and (49),

" # " #
0.9995 0 1.853e−10 0
Φ= and Qzd = .
0 1 0 1.960e−10

20 3) Using (22), the continuous-time measurement noise PSD is SN = 1.089e−5 m2 /s3 . Using
(62), the discrete-time measurement noise covariance is Qηd = 1.089e−3 m2 /s4 .
22 4) Starting from an initial condition of zero, the state vector ⃗xz (k) is propagated using (44)
and the IMU error z(k) is computed using (45). The data frequency is 100 Hz, so that
24 T = 0.01 s. For simulation purpose, L = 107 data samples (or 105 s) were generated.

19
 5) Figure 9 shows the ASD plot for the simulated data as the red dashed curve. The blue x’s
2 show the ASD curve for the IMU data, which is the same as is shown in Figures 2 and 8.
The green curve is the ASD for the analytic model computed by (42), which is the same
4 as is shown in Figure 8.

The close match between the green and red curves verifies the continuous-time to discrete-time
6 state-space model transformation. The closeness of the fit between those two curves and the
actual ASD plot depends on the choice of the structure and parameters of the continuous-time
8 state-space model.

Discussion of Issues and Tradeoffs

10 The previous sections presented an example approach to extract a discrete-time state-space

model to approximate the IMU error characteristics, as quantified by an ASD graph. There is
12 no single correct method. This section discusses various related issues and tradeoffs.

User-Acquired Data

14 When the manufacturer does not provide an ASD plot (or provides information that
the designer considers insufficient), the designer may contemplate acquiring their own data
16 to construct the ASD plots. This data acquisition process should be carefully designed after
considering the appropriate technical specifications [28], [29].

18 The sensor model in (3) shows that the sensor reading is a function of deterministic
parameters, stochastic errors, and the actual signal u(t). The ASD plot is only intended to
20 characterize the stochastic errors z(t). The data acquisition setup to acquire data to produce
the ASD of the instrument must isolate the IMU from the environment well enough that the
22 contributions to the measurement ũ from u are negligible. Under the best circumstances, this
is done on the lowest subterranean floor of a building, with the instrument attached to a heavy
24 mass on a vibration isolation system. Less than optimal results can be expected for an instrument
placed on a desk in an office on a higher floor, due to motion of the building or vibrations in
26 the floor, for example.

Optimization-based Parameter Selection for a Given Model

28 Once the designer selects a model structure, it is possible to define an optimization problem
to select the model parameters. See Section 12.11.4.1.2 in [28].

30 For example, given the state-space model structure described in the “Cumulative Error

20
 Model: N , B, K” section, the ASD model corresponding to (42) is
SB TB2
  S τ
2 SN TB  − Tτ − T2τ K
σz (τ ) = + 1− 3 − 4e B + e B + ,
τ τ 2τ 3
2 which can be written using the (positive) parameter vector θ⃗ = [ SN SB SK TB ] as
2
  θ τ
2 ⃗ = θ 1 θ2 θ4 θ 4

− τ
− 2τ
3
σz (τ ; θ) + 1− 3 − 4e θ4 + e θ4 + . (63)
τ τ 2τ 3
A cost function such as
L
X
2
⃗ =
C(θ) wi σ̂u2 (τi ) − σz2 (τi ; θ) (64)
i=1

can be optimized over positive values of θ. ⃗ The cost C(θ)

⃗ is computed using the known values
of the σ̂u2 (τi ) for i = 1, . . . , L/2 as plotted in the ASD plot and the weights
1
wi = ,
σ 2 [σ̂u2 (τi )]
where σ 2 [σ̂u2 (τi )] is the variance of the computed σ̂u2 (τi ), approximated as [51]

σ 2 σ̂u2 (τi ) = (2 σ [σ̂u (τi )])2 ,

 

4
⃗ in (63) is linear in three parameters
with σ [σ̂u (τi )] defined in (15). The functional form of σz2 (τ ; θ)
and only nonlinear in one. Therefore, the nonlinear search is only over θ4 . For each value of θ4 ,
6 the other optimal parameter values can be explicitly computed.

The results of performing this optimization are shown in the fourth row of Table 2. Note
8 that the optimization-based approach decreased the size of B, which can be interpreted as the
bias instability errors being less important than thought by the designer who performed the
10 manual tuning. This might motivate the designer to study whether the two-state N , B, K model
actually provides performance improvement under application conditions relative to the one-state
12 N , K model. The engineering art is to select the appropriate model structure that enables the
state estimator to calibrate the IMU sufficiently well to achieve a specified performance given
14 the implementation tradeoffs for a particular application.

A few words of caution are appropriate regarding the interpretation of such optimization-
16 based approaches to parameter selection. Note that “optimization-based” is not necessarily
“optimal”. This approach yields the parameter vector that minimizes a cost function for a given
18 model structure and choice of weights. Different cost functions, weights, or a different model
will yield different “optimized” model parameters. Also, many designers view the IMU in their
20 application as not quite as good as the IMU that the manufacturer specified. Therefore, their
preference is to select the IMU stochastic error model to slightly overbound the specified AV
22 plot rather than to optimally match that curve.

21
 Observability: Where is the Bias?

2 It is often the case that one of the deterministic errors accounted for in (3) is a constant
unknown bias b. For clarity of the discussion in this section, if all other deterministic errors are
4 assumed to be zero, then (3) has the form

ũ(t) = u(t) + b + z(t). (65)

The detrending process [mentioned before (12) in reference to the computation of the AV],
6 estimates and removes a bias from the specific set of data used to create the ASD plot. Over
a set of experiments, the variance of this “turn-on bias” could be determined to have the value
8 Pb0 . Because the IMU turn-on bias may change from one run to the next, the unknown portion
is considered as a constant bias b with the differential equation

ḃ(t) = 0, (66)

10 with the initial covariance Pb (0) = Pb0 [3], [9], [56].

Nonetheless, the ASD plot may still dictate the inclusion of the state zK with model

żK (t) = ωK (t), (67)

12 with the PSD for ωK (t) being SK = K 2 and initial covariance PzK (0) = cov(zk (0)) = 0.

Including both b and zK in the augmented state vector is problematic, as can be concluded
from a simple observability analysis. Assume for a moment the ideal situation where the signal
y(t) = ũ(t) − u(t) is available. Also, assume that z(t) = zK (t) (accounting for additional
stochastic errors would not improve the observability of b and zK ), then (65)-(67) are equivalent
to
" # " #" #
ḃ(t) 0 0 b(t)
= (68)
żK (t) 0 0 zK (t)
" #
h i b(t)
y(t) = 1 1 , (69)
zK (t)
where the noise ωK (t) has been dropped because it is not relevant to an observability analysis.
For this system, the observability matrix is
" #
1 1
O= , (70)
0 0
which has rank equal to one. Further analysis shows that it is only possible to estimate the sum
14 b(t) + zK (t) (see Example 3.19 in [3]). Alternative modeling approaches are discussed in the
“Tradeoffs in State Augmentation” section.

22
 Computational Impact of Augmented States

2 The augmented state vector, as defined in (5), has dimension nx = nv + nd + nz . The

following discussion considers the dimension of ⃗xv ∈ ℜnv as fixed and known, typically nv = 9
4 [see (4)].

For IMU state augmentation, the designer must choose the definition and model structure
6 for the deterministic and stochastic IMU error state vectors ⃗xd ∈ ℜnd and ⃗xz ∈ ℜnz . Because
there are three gyros and three accelerometers in the IMU, if ng states are augmented per gyro
8 and na states per accelerometer, then (nd + nz ) = 3(na + ng ). Therefore, nx = nv + 3(na + ng ).
Increasing the number of augmented states allows for a higher fidelity model, but it comes
10 with tradeoffs. The computational load of the extended Kalman filter increases in proportion to
n3x [that is O(n3x )] see [5], [9], [56]. Also, the desired increase in performance expected to be
12 gained with the additional computations may not be realized. For example, as the number of
augmented states is increased, observability issues may arise. While it is clear that unobservable
14 states increase the computational load without providing any performance improvement, it may
be less clear that weakly observable states may cause the implementation to be less robust to
16 unmodeled dynamics and nonlinearities.

Tradeoffs in State Augmentation

18 In sensor fusion applications, the IMU measurements are not processed in the state
estimation (for example the Kalman filter) measurement update. Instead, they are treated as
20 known inputs and considered during the time propagation of the state vector [see (2)].

To illustrate the tradeoffs related to state augmentation for modeling stochastic IMU
22 errors in such a filtering context, the differential equation for velocity shall be considered. In
a navigation frame mechanization, this velocity differential equation is given by (see Section
24 11.24 in [3])
⃗v˙ n = Cn f⃗ b − (2 ω
eb b ib ⃗ n+ω
ie⃗ n ) × ⃗v n + ⃗g n .
en eb l (71)

In this notation, ⃗vebn ∈ ℜ3 is the velocity of the platform (that is, body frame) with respect to Earth
26 and represented in navigation frame coordinates, Cnb is the direction cosine matrix describing the
rotation from body frame to navigation frame, ω ⃗ ien is the Earth rotation rate, ω n
⃗ en is the platform
n
28 transport rate, and ⃗gl is the local gravity vector. The axes of the navigation frame point in the
directions north, east, and down. The axes of the body frame usually coincide with the sensitive
30 axes of the IMU. The quantity f⃗ibb is the specific force vector, which is defined as the difference
between platform acceleration and local gravity vectors. The accelerometer triad contained in
˜
32 the IMU provides measurements f⃗ibb of this specific force, which – assuming no deterministic

23
 errors – relates to the true specific force via the stochastic error vector ⃗za ∈ ℜ3 :
˜
f⃗ibb = f⃗ibb + ⃗za . (72)

2 Each scalar component of the vector ⃗za is modeled as in (6)-(7). Substituting (72) into (71)
yields
˜
⃗v˙ ebn = Cnb f⃗ibb − (2 ω
⃗ ien + ω n
⃗ en ) × ⃗vebn + ⃗gln − Cnb ⃗za . (73)

4 Consider two accelerometer error modeling scenarios:

1) If only velocity random walk errors would be considered (see the “Random Walk Errors:
6 Angle and Velocity: zN (t)” section), the Kalman filter state vector would not need to be
augmented. The velocity random walk simply ends up as process noise in the Kalman
8 filter system model without the need for additional states. This model could be accurate
for small time durations (τ smaller than a few seconds for the IMU in Figure 2), but
10 would largely ignore the other sensor characteristics described by the corresponding ASD
plot for larger τ . It would not allow the sensor fusion algorithm to calibrate the IMU (for
12 example, estimate and remove a bias) on the fly.
2) Using the first-order Gauss-Markov error model (or the acceleration random walk model)
14 would add one accelerometer “bias state” per axis. In this case, the IMU stochastic error
model can match the ASD plot out to several tens of seconds. The estimated bias state
16 allows removal of the effect of the sensor bias (that is, calibration).

For the second case, both options (that is, the first-order Gauss-Markov error model or the
18 acceleration random walk model) have their pros and cons. With an acceleration random walk
model, the variance of the bias state grows linearly with time (that is, without bound) during
20 time intervals when the bias is not observable (see Section 4.6.3.2 in [3]). The physical inertial
sensor bias is of course bounded. When the bias becomes observable through aiding measurement
22 information and vehicle motion, the unrealistic growth of the bias variance may cause the
estimator gain to be unreasonably large. With the Gauss-Markov model, the variance of the
24 bias state stays bounded, even during time intervals when the bias state is not observable (see
(4.102) in [3]). However, when the bias is unobservable, the bias state estimate itself tends to
26 zero during time-propagation (that is, the previously estimated quantity is slowly forgotten).
Whether this is relevant or not depends on the duration of time without observability of the bias
28 state and Gauss-Markov model correlation time.

With the Gauss-Markov choice, (20) reduces to

⃗za = ⃗zN + ⃗zG .

24
 The augmented system model is given by
"
⃗v˙ ebn
# " ˜ #
Cnb f⃗ibb − (2 ω
⃗ ien + ω n
⃗ en ) × ⃗vebn + ⃗gln − Cnb ⃗zG
=
⃗z˙G −⃗µB ⃗zG
" #" #
−Cnb ω
⃗N
+ , (74)
I ω
⃗B

2 with the velocity random walk (that is, ω ⃗ N ) and the driving noise of the Gauss-Markov model
(that is, ω
⃗ B ) constituting the process noise.

4 This state augmentation allows for a calibration of the inertial sensors during the mission,
which improves the inertial navigation performance, especially when no aiding information is
6 available. The time-correlated nature of the inertial sensor biases is respected, and the sensor
characteristics represented in the corresponding ASD plot are modeled to the extent (that is,
8 cluster duration) where they are meaningful for the application. For example, typically in a
GNSS/INS system, GNSS measurements are processed with a rate of approximately 1 Hz. This
10 means that the time interval over which the INS propagates the state, without aiding corrections,
is normally limited to 1 s. Using inertial sensors up to tactical grade, an inertial-only navigation
12 during GNSS outages is only meaningful for a few minutes at best, because the inertial navigation
errors grow with time relatively quickly. In consequence, from the perspective of the filter design,
14 it is not required to model inertial sensor characteristics that become dominant in the ASD plot
for time intervals of several hundreds or thousands of seconds.

16 While the advantages of such an augmentation have just been outlined, as discussed earlier,
the augmentation with additional states increases the computational load. Before augmentation
18 with inertial sensor biases, a tightly coupled GNSS/INS filter would typically require 11 states:
three position error states, three velocity error states, three attitude error states, and two states
20 for the receiver clock error model. For multi-GNSS implementations, additional clock model
states are required per constellation. The augmentation with accelerometer and gyroscope bias
22 states (either ⃗zG or ⃗zK ) adds six states, which increases the computational load by roughly a
factor of four. Augmentation with two states per instrument (for example, both ⃗zG and ⃗zK ) adds
24 12 states, which increases the computational load by roughly a factor of 10.

Deterministic Errors

26 So far (for stochastic error modeling), only the ASD plots provided by the sensor
manufacturer have been considered. However, when designing a sensor fusion filter, additional
28 “deterministic” aspects of the IMU may need to be addressed. Inertial sensor manufacturers
usually provide a variety of specifications in addition to the ASD plots, which are addressed

25
 briefly below.

2 Misalignment. The IMU axis-to-axis misalignment describes the non-orthogonality of the

sensitive axes of the sensors. Some manufacturers provide a typical standard deviation for the
4 angle by which two sensor axes may deviate from an ideal 90 deg. With sufficient dynamics
and inertial sensor accuracy, misalignment could be estimated in the navigation filter, adding
6 nine states defined as follows. The IMU sensitive axes span the body frame. When considering
misalignment, it has to be defined how this is to be understood. One possible approach is to
8 define the sensitive axis of the z-accelerometer as z-axis of the body frame. The nonorthogonality
of the y-accelerometer axis to the z-accelerometer axis can be described with a single angle.
10 The x-axis is defined as the unique vector that is orthogonal to the y-z plane, so that the body
frame is defined completely. The description of the misalignment of x-accelerometer, and x, y,
12 and z-gyroscopes with respect to this body frame requires two angles for each sensor. For three
sensors, this requires up to nine states total. An IMU axis to platform frame misalignment (which
14 is to be understood as a simple rotation of the body frame with respect to the IMU frame that
it is supposed to be aligned with) has no negative impact on inertial navigation performance.

16 Nonlinearity. This term refers to the deviation of the sensor input-output curve from a straight
line. A line is fit to the input-output curve and the maximum deviation of the input-output curve
18 with respect to this fitted line, divided by the sensor measurement range (full scale), is defined
as nonlinearity. The nonlinearity usually cannot be estimated in the navigation filter.

20 Scale Factor. In the definition of nonlinearity, the line would ideally have a slope of one. The
deviation of the actual slope from one is the scale factor error (often stated in parts per million).
22 Scale factor error could potentially be estimated in the navigation filter, adding six states: three
for the gyroscopes and three for the accelerometers. However, as for the misalignment, the
24 linear scale factor error is often difficult to observe. Without vehicle motion, both scale factor
and misalignment are indistinguishable from biases.

26 Linear acceleration effect. MEMS gyroscopes show an acceleration-dependent bias error. Usually,
a navigation filter does not consider this type of error explicitly.

28 Vibration rectification error. This is especially relevant to accelerometers that exhibit a vibration
dependent bias. Possible reasons for such a vibration dependent bias include nonlinearity and
30 aliasing. Usually, a navigation filter does not consider this type of error explicitly.

The above brief discussion of deterministic sensor errors states that some of them (for
32 example, nonlinearity) cannot be accounted for explicitly in the system model of the state
estimator. Whether the augmentation of states for scale factor and misaligment calibration is

26
 worthwhile is application dependent. First, not all of them are observable individually. Second,
2 the increase in computational complexity could be prohibitive. Furthermore, the IMU is often
exposed to vibrations. Vibrations are micro-movements of the sensor, which may not be resolved
4 correctly by the IMU and INS strapdown algorithm. For example, the frequency of some of the
vibrations may exceed the Nyquist frequency of the sensor, causing aliasing to occur. These
6 errors might be modeled as additional time-correlated noise. However as the vibrational model
for each sensor would require at least a second-order Gauss-Markov model, such an approach
8 would require twelve additional states. Instead of explicitly considering vibration induced noise,
misalignment, linear and nonlinear scale factor errors, an often used approach is to increase the
10 navigation filter process noise beyond the levels that have been obtained from an analysis of the
ASD plot.

12 IMU Manufacturer Terminology

This section briefly discusses various additional terms that appear on some manufacturer
14 data sheets.

Bias in-run stability. The bias in-run stability is the component of the total sensor bias that
16 varies with time in a correlated fashion. This correlated temporal variation is typically modeled
as a first-order Gauss-Markov process, that is, zG (k). The bias in-run stability (or in-run bias)
18 is described on IMU data sheets with units that correspond to the ASD. In accordance with the
IEEE specifications [28], [30], this value corresponds to B at the minimum of the ASD plot,
20 as discussed relative to (32). However, the cluster time corresponding to this minimum value
is typically not provided by sensor manufacturers. It can be extracted from the ASD plot, if
22 provided; otherwise, the correlation time must be selected and tuned by the designer.

Turn-on bias. The turn-on bias (also referred to as the start-up, run-to-run, or repeatability bias)
24 is an offset in the sensor readings, that potentially changes each time the sensor is switched on.
The total inertial sensor bias can be seen as the sum of this turn-on bias and a time varying
26 contribution (that is, bias in-run stability). The turn-on bias is relevant for the initialization of the
variance of the inertial sensor bias states in the navigation filter. Using the state vector defined
28 relative to (48) as an example, the initial covariance of the augmented states zG (0) and zK (0)
should be selected to add up to the manufacturer specified variance of the turn-on bias. As
30 explained in the “Observability: Where is the Bias?” section, the turn-on bias cannot be inferred
from an ASD plot.

Angle/velocity random walk. The random walk describes the impact of the sensor inherent
white noise when integrating. The unit is √◦h for gyroscopes and m/s
√
h
for accelerometers. An

27
 angular random walk of α √◦h indicates that after integrating angular rate measurements for H
√
hours, an angle error standard deviation of α H degrees will result due to the sensor inherent
white noise. The random walk contributes to the process noise in the navigation filter system
model [for example, see (74)]. The angle/velocity random walk parameter is the square root of
the PSD of the sensor inherent white noise, which is denoted by N throughout this article. In
some data sheets, the angle/velocity random walk is referred to as the noise density, typically
m/s2
provided with the unit √◦/hHz
for gyroscopes and √
Hz
for accelerometers. The noise density is
an alternative representation to the random walk parameter N . It is possible to convert between
both representations of the sensor inherent white noise. For example,
√
m/s2 m s m/s m/s
√ = = √ = 60 √ .
Hz s2 s h

Available Software Package

2 An open-source Matlab software package is available at https://github.com/jaffarrell/AV-

Matlab-SW. The software has two components. The first component completes the following: (1)
4 Given a continuous-time error model, it computes an equivalent discrete-time state-space model.
(2) Given that discrete-time model, a simulation produces a stochastic error sequence zk = z(tk ),
6 where tk = k T for k = 1, . . . , L. (3) Given a sequence of stochastic errors {zk }Lk=1 , it computes
p
the AV (i.e., σu2 (τ )) and plots the ASD (i.e. σu2 (τ )). The second component implements the
8 optimization-based approach described in (63)-(64) to fit the parameters of the NBK continuous-
time state-space model described in (40)-(42) to a given set of AV data. Together, these two
10 components enable a complete design cycle. Inclusion of the optimization-based approach is for
completeness; it is not meant to imply that it is the recommended approach for selecting the
12 state-space stochastic error model.

This software release includes a dataset for the demonstration example herein. The user
14 must adapt that approach and software to the model appropriate for their instrument.

Conclusions

16 The main purpose of this article has been to present a tutorial on the process that starts
from an instrument’s ASD plot (or its derived parameters) and constructs a state-space IMU
18 error model suitable for real-time INS error state estimation and IMU calibration using data
fusion methods such as the KF, EKF, UKF, PF, or MAP. An example model construction and
20 verification method is included. We do not claim that this approach is unique or optimal. It is
representative of industry-standard methods.

28
 In addition, this tutorial included an extensive discussion of the issues and tradeoffs that a
2 designer should consider, including performance, computational load, observability, and extent
of the cluster time τ that is relevant to a given application.

4 Finally, the manufacturer provided ASD plot or parameters should be considered as a

starting point. It will dictate the dominant forms of error and reasonable values for the error model
6 parameters. Minor tuning relative to those reasonable parameters will normally be required to
accommodate the particular IMU that is available, and error terms that were ignored or neglected
8 due to observability or computational reasons. A considerable tuning might be required, in case
the IMU is exposed to vibrations containing frequencies exceeding the Nyquist frequency.

10 Acknowledgments

The authors greatly appreciate the efforts of the anonymous reviewers and editorial team.
12 Their comments enhanced the presentation and readability of the article. Specifically, the model
robustness ideas of the penultimate sentence in the penultimate paragraph of the “Bias Instability:
14 zB (t)” section and in the last sentence in the “Computational Impact of Augmented States”
section are based on suggestions from the reviewers.

16 The authors also thank Zhiqiang Xing and Demoz Gebre-Egziabher for supplying the IMU
data and answering questions that helped to create Figure 1; and Y. Yang of Nav Technology
18 for supplying the IMU data that is the basis of Figures 2, 8, and 9.

The optimization-based approach is not original to this article. It is practiced in industry.

20 The general idea was explained to the first author by John R. Dowdle and Karl W. Flueckiger
in the early 1990’s.

22 This work was supported in part by the Coordination for the Improvement of Higher
Education Personnel (CAPES), under grant 88881.169927/2018-01; the Brazilian Agricultural
24 Research Corporation (EMBRAPA), under grant 212-20/2018; and, the Brazilian National
Council for Scientific and Technological Development (CNPq), under grant 313160/2019-8.

29
 References

2 [1] K. R. Britting, Inertial Navigation Systems Analysis. John Wiley & Sons Canada, Limited,
1971.
4 [2] M. M. Kuritsky and M. S. Goldstein, “Inertial navigation,” Proc. of the IEEE, vol. 71,
no. 10, p. 1156–1176, 1983.
6 [3] J. A. Farrell, Aided Navigation: GPS with High Rate Sensors. McGraw-Hill, Inc., 2008.
[4] P. D. Groves, Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems.
8 London: Artech House Remote Sensing Library, 2013.
[5] M. S. Grewal, A. P. Andrews, and C. G. Bartone, Global Navigation Satellite Systems,
10 Inertial Navigation & Integration. John Wiley & Sons, 2013.
[6] D. A. Grejner-Brzezinska, C. K. Toth, T. Moore, J. F. Raquet, M. M. Miller, and A. Kealy,
12 “Multisensor navigation systems: A remedy for GNSS vulnerabilities?” Proc. of the IEEE,
vol. 104, no. 6, p. 1339–1353, 2016.
14 [7] R. E. Kalman, “A new approach to linear filtering and prediction problems,” J. Basic Eng.
- Trans. ASME, vol. 82, no. D, pp. 35–45, 1960.
16 [8] R. E. Kalman and R. S. Bucy, “New results in linear filtering and prediction theory,” J.
Basic Eng. - Trans. ASME, vol. 1, pp. 95–108, 1961.
18 [9] A. Gelb, J. F. Kasper Jr, R. A. Nash Jr, C. F. Price, and A. A. Sutherland Jr, Applied
Optimal Estimation. Cambridge: Analytic Sciences Coorporation, 1974.
20 [10] F. M. Mirzaei and S. I. Roumeliotis, “A Kalman filter-based algorithm for IMU-camera
calibration: Observability analysis and performance evaluation,” IEEE Trans. Robot.,
22 vol. 24, no. 5, pp. 1143–1156, 2008.
[11] R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman
24 Filtering with Matlab Exercises, 2nd ed. Wiley, 1992.
[12] J. A. Farrell, T. Givargis, and M. J. Barth, “Real-time differential carrier phase GPS-aided
26 INS,” IEEE Transactions on Control Systems Technology, vol. 8, no. 4, pp. 709–721, 2000.
[13] A. I. Mourikis and S. I. Roumeliotis, “A Multi-State Constraint Kalman Filter for Vision-
28 aided Inertial Navigation,” IEEE Int. Conf. Robotics and Automation, 2007.
[14] J. Civera, A. J. Davidson, and J. M. Montiel, “Inverse depth parametrization for monocular
30 SLAM,” IEEE T. Robotics, vol. 24, no. 5, pp. 932–945, 2008.
[15] M. Li, B. H. Kim, and A. I. Mourikis, “Real-time motion tracking on a cellphone using
32 inertial sensing and a rolling-shutter camera,” IEEE Int. Conf. Robotics and Automation,
2013.
34 [16] S. Julier, J. Uhlmann, and H. Durrant-Whyte, “A new method for the nonlinear trans-
formation of means and covariances in filters and estimators,” IEEE T. Autom. Control,
36 vol. 45, no. 3, p. pages, 477–482.

30
 [17] S. Julier and J. Uhlmann, “Unscented filtering and nonlinear estimation,” Proc. of the
2 IEEE, vol. 92, no. 3, p. 401–422, 2004.
[18] J. L. Crassidis, “Sigma-point Kalman filtering for integrated GPS and inertial navigation,”
4 IEEE T. Aerospace and Electronic Systems, vol. 42, no. 2, pp. 750–756, 2006.
[19] J. S. Liu and R. Chen, “Sequential Monte Carlo methods for dynamic systems,” J. of the
6 American Stat. Assoc., vol. 93, no. 443, pp. 1032–1044, 1998.
[20] M. Pitt and N. Shephard, “Filtering via simulation: auxiliary particle filter,” J. of the
8 American Stat. Assoc., vol. 94, pp. 590–599, 1999.
[21] A. Doucet, J. de Freitas, and N. Gordon, Eds., Sequential Monte Carlo Methods In
10 Practice. Springer-Verlag, 2001.
[22] S. Thrun, W. Burgard, and D. Fox, “A probabilistic approach to concurrent mapping
12 and localization for mobile robots,” Mach. Learn., vol. 31, no. 1-3, pp. 29–53, 1998.
[Online]. Available: http://dx.doi.org/10.1023/A:1007436523611
14 [23] F. Dellaert and M. Kaess, “Square Root SAM: Simultaneous localization and mapping via
square root information smoothing,” Intl. J. of Robotics Research, IJRR, vol. 25, no. 12,
16 pp. 1181–1204, Dec 2006.
[24] R. Eustice, H. Singh, J. J. Leonard, M. Walter, and R. Ballard, “Visually navigating the
18 RMS Titanic with SLAM information filters,” Robotics: Science and Systems, 2005.
[25] R. Eustice, H. Singh, and J. J. Leonard, “Exactly sparse delayed-state filters for view-based
20 SLAM,” IEEE T.-Rob., vol. 22, p. 1100 – 1114, 2006.
[26] M. Kaess, A. Ranganathan, and F. Dellaert, “iSAM: Incremental Smoothing and Mapping,”
22 IEEE T.-Rob., vol. 24, no. 6, p. 1365 – 1378, 2008.
[27] S. Zhao, Y. Chen, H. Zhang, and J. A. Farrell, “Differential GPS aided inertial navigation:
24 A contemplative realtime approach,” IFAC, pp. 8959–8964, 2014.
[28] IEEE, “Specification Format Guide and Test Procedure for Single-Axis Interferometric
26 Fiber Optic Gyros,” Standard 952, R2003.
[29] ——, “Specification Format Guide and Test Procedure for Linear, Single-Axis, Non-
28 Gyroscopic Accelerometers,” Standard 1293, 1999.
[30] ——, “Specification Format Guide and Test Procedure for Single-Axis Laser Gyros,”
30 Standard 647, R2006.
[31] J. Baziw and C. T. Leondes, “In-flight alignment and calibration of inertial measurement
32 units - Part I: General formulation,” IEEE Trans. Aero. Elec. Sys., vol. 8, no. 4, pp.
439–449, 1972.
34 [32] A. J. Dierendonck, J. B. McGraw, and R. G. Brown, “Relationship between Allan variances
and Kalman filter parameters,” in Proc. 16th Annual Precise Time and Time Interval (PTTI)
36 Applications and Planning Meeting. Greenbelt, MD: NASA Goddard Space Flight Center,
1984, pp. 273–293.

31
 [33] J. J. Ford and M. E. Evans, “Online estimation of Allan variance parameters,” Journal of
2 Guidance, Control, and Dynamics, vol. 23, no. 6, pp. 980–987, 2000.
[34] N. El-Sheimy, H. Hou, and X. Niu, “Analysis and modeling of inertial sensors using
4 Allan variance,” IEEE T. on Instrumentation and Measurement, vol. 57, no. 1, pp. 140–
149, 2008.
6 [35] X. Zhang, Y. Li, P. Mumford, and C. Rizos, “Allan variance analysis on error characters
of mems inertial sensors for an FPGA-based GPS/INS system,” International Symposium
8 on GPS/GNSS, p. 127–133, 2008.
[36] Z. Xing and D. Gebre-Egziabher, “Modeling and bounding low cost inertial sensor errors,”
10 in Proc. Position Location and Navigation Symp. Monterey: IEEE, 2008, pp. 1122–1132.
[37] V. Saini, S. C. Rana, and M. M. Kube, “Online estimation of state space error model for
12 MEMS IMU,” Journal of Modelling & Simulation of Systems, vol. 1, no. 4, pp. 219–225,
2010.
14 [38] V. Vaccaro and A. Zaki, “Statistical modeling of rate gyros,” IEEE T. Instrumentation and
Measurement, vol. 61, no. 3, pp. 673 – 684, 2012.
16 [39] J. Hidalgo-Carrió, S. Arnold, and P. Poulakis, “On the design of attitude-heading reference
systems using the Allan variance,” IEEE T. Ultrason., Ferroelectr., Freq. Control, vol. 63,
18 no. 4, pp. 656–665, 2016.
[40] F. O. Silva, E. M. Hemerly, and W. C. Leite Filho, “On the error state selection for sta-
20 tionary SINS alignment and calibration Kalman filters - Part II: Observability/estimability
analysis,” Sensors, vol. 17, no. 439, pp. 1–34, 2017.
22 [41] D. H. Titterton and J. L. Weston, Strapdown Inertial Navigation Technology. Reston:
Institution of Electrical Engineers, 2004.
24 [42] R. Zhang, F. Hoflinger, and L. M. Reindl, “Calibration of an IMU using 3-D rotation
platform,” IEEE Sens. J., vol. 14, no. 6, pp. 1778–1787, 2014.
26 [43] Y. Stebler, S. Guerrier, J. Skaloud, and M. P. Victoria-Feser, “A framework for inertial
sensor calibration using complex stochastic error models,” in Proc. Position, Location and
28 Navigation Symp. Myrtle Beach: IEEE, 2012, pp. 849–861.
[44] A. G. Quinchia, G. Falco, E. Falletti, F. Dovis, and C. Ferrer, “A comparison between
30 different error modeling of MEMS applied to GPS/INS integrated systems,” Sensors,
vol. 13, no. 1, pp. 9549–9588, 2013.
32 [45] Crossbow, “Specification sheet for the µnav navigation & servo control board,” Document
No. 6020-0083-02, Rev. B., 2006.
34 [46] S. Nassar, K. P. Scharwz, and N. El-Sheimy, “Modeling inertial sensor errors using
autoregressive AR models,” Navigation, vol. 51, no. 4, pp. 259–268, 2004.
36 [47] Z. Miao, F. Shen, D. Xu, K. He, and C. Tian, “Online estimation of Allan variance

32
 coefficients based on a neural-extended Kalman filter,” Sensors, vol. 15, no. 1, pp. 2496–
2 2524, 2015.
[48] W. J. Riley, Handbook of Frequency Stability Analysis, ser. Volume 1065 of NIST
4 special publication. U.S. Department of Commerce, National Institute of Standards and
Technology, 2008.
6 [49] D. W. Allan, “Statistics of atomic frequency standards,” Proc. of the IEEE, vol. 54, no. 2,
pp. 221–230, 1966.
8 [50] J. A. Barnes, A. R. Chi, L. S. Cutler, D. J. Healey, D. B. Leeson, T. E. McGunigal,
J. A. Mullen Jr., W. L. Smith, R. L. Sydnor, R. F. C. Vessot, and G. R. Winkler,
10 “Characterization of frequency stability,” IEEE Trans. Instrum. Meas., vol. IM-20, no. 2,
pp. 105–120, 1971.
12 [51] P. Lesage and C. Audoin, “Characterization of frequency stability: uncertainty due to the
finite number of measurements,” IEEE Trans. Instrum. Meas., vol. 22, no. 5, pp. 157–161,
14 1973.
[52] ——, “Correction to ”characterization of frequency stability: uncertainty due to the finite
16 number of measurements”,” IEEE Trans. Instrum. Meas., vol. 25, no. 3, pp. 270–270,
1976.
18 [53] M. M. Tehrani, “Ring laser gyro data analysis with cluster sampling technique,” Proc.
SPIE, vol. 412, no. 1, pp. 207–220, 1983.
20 [54] J. Li and J. Fang, “Not fully overlapping Allan variance and total variance for inertial
sensor stochastic error analysis,” IEEE Trans. Instrum. Meas., vol. 62, no. 10, pp. 2659–
22 2672, 2013.
[55] H. Hou, Modeling Inertial Sensors Errors Using Allan Variance. MS Thesis, University
24 of Calgary, September 2004, no. UCGE 20201.
[56] R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman
26 Filtering with Matlab Exercises. Wiley, 2012.
[57] W. Gautschi and W. F. Cahill, Exponential Integral and Related Functions: Handbook of
28 Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York:
Dover, 1972.
30 [58] P. G. Savage, Strapdown Analytics. Maple Plain: Strapdown Associates, 2007.
[59] N. J. Kasdin, “Discrete simulation of colored noise and stochastic processes and 1/f α
32 power law noise generation,” Proc. of the IEEE, vol. 83, no. 5, pp. 802–827, 1995.
[60] C. Van Loan, “Computing integrals involving the matrix exponential,” IEEE Transactions
34 on Automatic Control, vol. 23, no. 3, 1978.
[61] W. Wrigley, “History of inertial navigation,” NAVIGATION: Journal of the Institute of
36 Navigation, vol. 24, no. 1, Spring 1977.
[62] A. Lawrence, Modern Inertial Technology. Springer, 1993.

33
 [63] P. Teunissen and O. Montenbruck, Eds., Springer Handbook of Global Navigation Satellite
2 Systems. Springer, 2017.
[64] P. S. Maybeck, Stochastic Models, Estimation and Control. London: Academic Press,
4 Inc., 1979.
[65] F. M. Ham and R. G. Brown, “Observability, eigenvalues, and Kalman filtering,” IEEE
6 Trans. Aero. Elec. Sys., vol. 19, no. 2, pp. 269–273, 1983.
[66] J. E. Potter and R. G. Stern, “Statistical filtering of spacecraft navigation measurements,”
8 Proc. AIAA GUidance and Control Conf., 1963.
[67] G. H. Golub, “Numerical methods for solving least squares problems,” Numer. Math.,
10 vol. 7, pp. 206–216, 1965.
[68] S. F. Schmidt, “Computational techniques in kalman filtering,” NATO-AGARD-139, 1970.
12 [69] ——, “Sequential square root filtering and smoothing of discrete linear systems,” Auto-
matica, vol. 10, pp. 147–158, 1974.
14 [70] K. R. Britting, “Self-alignment techniques for strapdown inertial navigation systems with
aircraft application,” J. Aircraft, vol. 7, no. 4, pp. 302–307, 1970.
16 [71] V. Krishnan and K. Grobert, “Initial alignment of a gimballess inertial navigation system,”
IEEE Trans. Automat. Contr., vol. 15, no. 6, pp. 667–671, 1970.
18 [72] J. Baziw and C. T. Leondes, “In-flight alignment and calibration of inertial measurement
units - part II: Experimental results,” IEEE Trans. Aero. Elec. Sys., vol. 8, no. 4, pp.
20 450–465, 1972.
[73] I. Y. Bar-Itzhack and B. Porat, “Azimuth observability enhancement during inertial
22 navigation system in-flight alignment,” J. Guid. Control, vol. 3, no. 4, pp. 337–344, 1980.
[74] B. Porat and B.-I. I. Y., “Effect of acceleration switching during INS in-flight alignment,”
24 J. Guid. Control, vol. 4, no. 4, pp. 385–389, 1981.
[75] I. Y. Bar-Itzhack and N. Berman, “Control theoretic approach to inertial navigation
26 systems,” J. Guid. Control, vol. 11, no. 3, pp. 237–245, 1988.
[76] D. Goshen-Meskin and I. Y. Bar-Itzhack, “Observability studies of inertial navigation
28 systems,” J. Guid. Control, vol. 11, no. 3, pp. 1283–1289, 1989.
[77] S. Hong, M. H. Lee, H. H. Chun, S. H. Kwon, and J. L. Speyer, “Observability of error
30 states in GPS/INS integration,” IEEE Trans. Aero. Elec. Sys., vol. 54, no. 2, pp. 731–743,
2005.
32 [78] I. Rhee, M. F. Abdel-Hafez, and J. L. Speyer, “Observability of an integrated GPS/INS
during maneuvers,” IEEE Trans. Aero. Elec. Sys., vol. 40, no. 2, pp. 1362–1366, 2004.
34 [79] F. O. Silva, E. M. Hemerly, and W. C. Leite Filho, “Error analysis of analytical coarse
alignment formulation for stationary SINS,” IEEE Trans. Aero. Elec. Sys., vol. 52, no. 4,
36 pp. 1777–1796, 2016.
[80] F. O. Silva, E. M. Hemerly, W. C. Leite Filho, and H. K. Kuga, “A fast in-field coarse

34
 alignment and bias estimation method for stationary intermediate-grade IMUs,” IEEE
2 Trans. Instrum. Meas., vol. 67, no. 4, pp. 831–838, 2018.
[81] A. E. Bryson, “Kalman filter divergence and aircraft motion estimators,” J. Guidance and
4 Control, vol. 1, no. 1, pp. 71–79, 1978.
[82] A. H. Jazwinski, “Limited memory optimal filtering,” IEEE TAC, vol. 13, no. 5, pp. 558–
6 563, 1968.
[83] J. H. L. K. R. Muske, J. B. Rawlings, “Receding horizon recursive state estimation,” ACC,
8 pp. 900–904, 1993.
[84] P. E. Moraal and J. W. Grizzle, “Observer design for nonlinear systems with discrete–time
10 measurements,” IEEE TAC, vol. 40, no. 3, p. 395–404, 1995.
[85] W. Li, W. Li, X. Cui, S. Zhao, and M. Lu, “A tightly coupled rtk/ins algorithm
12 with ambiguity resolution in the position domain for ground vehicles in harsh urban
environments,” Sensors, vol. 18, no. 7, p. 2160, 2018.
14 [86] E. Parzen, Stochastic Processes. San Francisco: Holden-Day, 1962.
[87] D. W. Allan and J. A. Barnes, “A modified ”Allan variance” with increased oscillator
16 characterization ability,” in Proc. 35th Annual Frequency control Symposium. Fort
Monmouth: USA ERADCON, 1981, pp. 471–475.
18 [88] D. A. Howe, D. W. Allan, and J. A. Barnes, “Properties of signal sources and measurement
methods,” in Proc. 35th Annual Frequency control Symposium. Fort Monmouth: USA
20 ERADCON, 1981, pp. A1–A47.
[89] J. Rutman, “Characterization of phase and frequency instabilities in precision frequency
22 sources: fifteen years of progress,” Proc. of the IEEE, vol. 66, no. 9, pp. 1048–1075, 1978.
[90] F. L. Walls and D. W. Allan, “Measurements of frequency stability,” Proc. of the IEEE,
24 vol. 74, no. 1, pp. 162–168, 1986.
[91] A. H. Jazwinski, Stochastic Processes and Filtering Theory, ser. Mathematics in Science
26 and Engineering. New York: Academic Press, Inc., 1970.
[92] R. E. Kalman, T. S. Englar, and R. S. Bucy, “Fundamental study of adaptive control
28 systems,” ASD-TR-61-27, vol. 1, 1962.
[93] P. Lesage and C. Audoin, “Estimation of the two-sample variance with a limited number
30 of data,” in Proc. 31st Annual Symposium on Frequency Control. Atlantic City: IEEE,
1977, pp. 311–318.
32 [94] K. Yoshimura, “Characterization of frequency stability: uncertainty due to the autocorre-
lation of the frequency fluctuations,” IEEE Trans. Instrum. Meas., vol. 27, no. 1, pp. 1–7,
34 1978.
[95] P. Lesage and C. Audoin, “Characterization of measurement of time and frequency
36 stability,” Radio Science, vol. 14, no. 4, pp. 521–539, 1979.
[96] C. A. Greenhall, “A structure function representation theorem with applications to

35
 frequency stability estimation,” IEEE Trans. Instrum. Meas., vol. 32, no. 2, pp. 364–370,
2 1983.
[97] D. Yuan, X. Ma, Y. Liu, Z. Shang, and S. Yan, “Statistical modeling of random walk
4 errors for triaxial rate gyros,” IEEE Trans. Instrum. Meas., vol. 65, no. 2, pp. 286–296,
2016.
6 [98] P. Lesage and C. Audoin, “Commetns on ”characterization of frequency stability: uncer-
tainty due to the finite number of measurements”,” IEEE Trans. Instrum. Meas., vol. 24,
8 no. 1, pp. 86–86, 1975.
[99] ——, “Effect of dead-time on the estimation of the two-sample variance,” IEEE Trans.
10 Instrum. Meas., vol. 28, no. 1, pp. 6–10, 1979.
[100] K. Yoshimura, “Degrees of freedom of the estimate of the two-sample variance in the
12 continuous sampling method,” IEEE Trans. Instrum. Meas., vol. 38, no. 6, pp. 1044–
1049, 1989.
14 [101] C. A. Greenhall, “Recipes for degrees of freedom of frequency stability estimators,” IEEE
Trans. Instrum. Meas., vol. 40, no. 6, pp. 994–999, 1991.
16 [102] T. Walter, “Characterizing frequency stability: A continuous power-law model with discrete
sampling,” IEEE Trans. Instrum. Meas., vol. 43, no. 1, pp. 69–79, 1994.
18 [103] P. Lesage and T. Ayi, “Characterization of frequency stability: Analysis of the modified
Allan variance and properties of its estimate,” IEEE Trans. Instrum. Meas., vol. 33, no. 4,
20 pp. 332–336, 1984.

36
 Sidebar: Nontechnical Article Summary

2 Autonomous vehicles utilize control systems to cause the vehicle state to follow a
desired trajectory. The control system incorporates information about the vehicle position,
4 velocity, acceleration, attitude, and angular rate, which can be computed by integration of the
measurements of an inertial measurement unit (IMU). This integrative process also accumulates
6 IMU measurement errors that cause the integrated IMU measurements (which is the vehicle
state estimate) to slowly diverge from the true vehicle state. The difference is referred to as the
8 vehicle error state vector. This vehicle error state vector and various IMU calibration parameters
can be estimated through the state estimation process using information from external sensors
10 [for example, camera, global navigation satellite systems (GNSS), Lidar, Radar]. Such real-time
calibration of the IMU (and other sensors) results in improved accuracy and slower rates of IMU
12 error accumulation during the time intervals between the measurements from external sensors.

Design of the state estimator requires definition of the IMU error state vector and
14 its stochastic discrete-time state-space model. This modeling process begins from the IMU
performance specification information provided by the IMU manufacturer, which (per the IEEE
16 standards) is communicated through the Allan Variance. Information extracted from the Allan
standard deviation graph allows the analyst to evaluate various issues and tradeoffs involved in
18 selecting the continuous-time IMU error state-space model. This tutorial article discusses this
approach, the issues and tradeoffs, the translation of the continuous-time model to an equivalent
20 discrete-time model for implementation, and a verification approach. Example Matlab scripts are
supplied that implement each step of this process.

37
 Sidebar: Aided Inertial Navigation System History

2 The first use of inertial navigation dates back to the German V-2 missile in 1942. After
World War II, the United States started to develop inertial navigation systems (INS) for ballistic
4 missiles, and (later in the 1960s) for the Apollo missions, as well as military and commercial
airplanes [61]. The early inertial navigation systems used a gimballed platform that decoupled
6 the host vehicle attitude changes from those of the platform. This decoupling allowed three
accelerometers mounted on the platform to maintain alignment with the north, east, and down
8 directions. For many applications, a reduced set of sensors was sufficient, for example, the vertical
axis and respective sensors were omitted. Because the initial position and velocity were known,
10 integration of the accelerometer measurements propagated the position and velocity vectors
forward in time. Gyroscopes were used to measure small perturbations of the platform attitude,
12 which were then compensated mechanically to maintain the alignment of the accelerometers
with the north, east and down directions.

14 These gimballed systems, despite their impressive accuracy, had several drawbacks. The
mechanical construction was complex, bulky, and expensive. Furthermore, aerobatic maneuvers
16 could cause a gimbal lock: When two gimbal axes become aligned, a rotation perpendicular to
this axis causes the platform to lose alignment. A solution was to add a fourth gimbal, which
18 further increased complexity and cost [62].

Advances in electronic and gyroscope technology in the 1960s enabled the development
20 of strapdown INS. In a strapdown INS, an inertial measurement unit (IMU) consisting of three
accelerometers and three gyroscopes is attached rigidly to the host vehicle. The gyroscopes
22 measure the attitude changes of the host vehicle, thereby enabling computational tracking of its
attitude. This requires that these gyroscopes be capable of accurately measuring angular rates of
24 up to several hundred degrees per second for highly agile host platforms, whereas for a gimballed
system, a gyroscope measurement range of few degrees per hour would have been sufficient.
26 Knowledge of the host vehicle attitude allows transformation of the accelerometer measurements
of the specific force from the body frame to an Earth reference frame (for example, a north-
28 east-down fixed tangent frame) in which compensation of gravity is straightforward; thereby
yielding measurements of the host vehicle acceleration. As in a gimballed system, this host
30 vehicle acceleration is then integrated twice to compute the changing host vehicle position and
velocity over time.

32 For any INS (gimballed or strapdown) the errors in the inertial measurements accumulate,
and together with the inaccuracies of the host vehicle initial state estimate, cause a deterioration
34 of the accuracy of the navigation solution with time. To counter this temporal growth of

38
 position, velocity and attitude errors, the inertial navigation system can be corrected using aiding
2 measurements from external sensors (for example, GNSS receiver, camera, lidar, radar).

A typical GNSS-aided INS architecture is shown in Figure 10. The strapdown algorithm
4 numerically integrates the nonlinear kinematic model described in (2) using the calibrated
IMU measurements ⃗uˆ(t) as inputs. This integration is at the sampling rate of the IMU (for
6 example, hundreds or thousands of times per second), which is designed to be high relative
to the bandwidth of both the IMU and the host vehicle. When aiding sensor measurements
8 are available, the navigation filter – often an extended or linearized Kalman filter – estimates
corrections for the state and calibration factors. Aiding measurements are typically available
10 at low rates (for example, one vector measurement per second), relative to the host vehicle
bandwidth. Calibration factors may include deterministic errors (for example, scale factor and
12 sensor-axis alignment errors) and time-correlated stochastic errors. This on-the-fly estimation of
IMU calibration factors continuously recalibrates the INS, leading to a superior INS performance
14 during time intervals when aiding sensor measurements are not available (for example, GNSS
outages).

16 Two aided INS system architectures are widely used: loose and tight coupling (see Ch. 28
in [63]). In a loosely coupled GNSS-aided approach, the GNSS receiver computes position and
18 velocity estimates internally from its pseudorange and Doppler measurements without using INS
information. This GNSS computation can occur only when the receiver has at least four satellites
20 in view. When the GNSS position and velocity measurements are output from the receiver, the
residual between them and their INS computed values are used to drive the navigation filter to
22 estimate the INS error state; otherwise, the INS continues to integrate without correction. In
a tightly coupled approach, the residual of the navigation filter is formed between the GNSS
24 pseudorange and Doppler measurements and the predictions of those quantities as computed
by the INS. This allows aiding even when fewer than four satellites are in view. The main
26 tradeoff is that tightly coupled systems offer the potential for higher-performance, especially
when fewer measurements are available, but are more complex to implement due to the need
28 for the navigation system: to process GNSS ephemeris data; to calculate satellite positions and
velocities; and, to apply corrections for ionosphere, troposphere, satellite clocks and broadcast
30 group delays (see Section 28.2 of [63]). Similar tradeoffs apply for aiding with alternative aiding
sensors.

32 There is a rich literature concerning dynamics, control, and system theoretic contributions
within the inertial navigation context [1], [64]. Throughout their history, the Kalman filter has
34 played an important role [7], [8], [65]. This history has included a focused effort on numerical
methods [66], [67], [68], [69]. Observability studies for both stationary and time-varying systems

39
 systems is critical both for initialization and on-the-fly calibration [31], [70], [71], [72], [73],
2 [74], [75], [76], [77], [78], [79], [80]. Also understanding the controllability of the state from
the perspective of the driving noise was critical to removing issues of Kalman filter divergence
4 [81]. More recently (but still with a long history [82], [83], [84]), efficient numeric methods for
real-time trajectory (as opposed to state) estimation (that is, real-time smoothing) are important
6 for inertial-based simultaneous location and mapping applications [23], [24], [25], [26].

The ongoing decreases in the cost of inertial sensors, aiding measurements, and compu-
8 tation are allowing aided strapdown INS to be feasible with respect to both cost and accuracy
in commercial applications. For example, the interest in combining GNSS-aided INS with Real-
10 Time Kinematic (RTK) techniques (see Chapter 26 in [63]) capable of achieving submeter
accuracy is growing. An overview on possible architectures is given in [85]. The performance
12 of aided inertial systems is dependent on the navigation filter incorporating a state-space model
of the IMU stoshastic errors. This article provides a tutorial describing the industry standard
14 process and tradeoffs related to defining such models.

40
 Sidebar: Simplified Inertial Navigation System Example

2 This section presents a simplified two-dimensional inertial navigation system (INS)

example. The purpose is to present a nonlinear kinematic model, the INS equations that would
4 propagate the vehicle state through time, the linearized model that predicts the growth in the
INS error over time, how that error model is used to propagate the error covariance through
6 time, and state augmentation for sensor calibration.

Various reference frames and variables are illustrated in Figure 11. Earth is assumed to be
8 circular with radius R, nonrotating, and uniform density. The point P represents the location of
the inertial measurement unit (IMU) on the vehicle. The vehicle is free to translate in ℜ2 and
10 rotate with one degree of freedom, denoted by θ. The IMU sensitive axes are in the directions
indicated by the unit vectors ⃗u and w, ⃗ which define the vehicle reference frame. The origin of
12 the geographic reference frame is defined as the projection of the point P onto the Earth surface
along the vector to the Earth’s center. The instantaneous Earth tangent plane at the origin of the
14 geographic frame defines the unit vectors ⃗n and d. ⃗ The height of the point P above the tangent
plane is the altitude h. The latitude ϕ and pitch θ are defined as positive in the directions indicated
16 in the figure. The vectors ⃗z and p⃗ define the axes of the Earth-centered reference frame.

The kinematic model for the IMU at point P (that is, the vehicle model) is
" # " #" #
1
ϕ̇ R+h
0 vn
= (75)
ḣ 0 −1 vd
" # " #
vn vd
v̇n g R+h
= ⃗aiv + −vn2 (76)
v̇d R+h

θ̇ = ωgv , (77)

18 where ⃗agiv is the inertially-referenced vehicle acceleration vector represented in the geographic
frame and ωgv is the rotation rate of the vehicle relative to the geographic frame. This angular
20 rate is computed as ωgv = ωiv − ωig , where ωiv is the rotation rate of the vehicle with respect to
vn
an inertial frame (which is measured by the gyro) and ωig = −ϕ̇ = − R+h is the transport rate
22 of the geographic frame with respect to the inertial frame. The superscript g (or v) on the vector
quantities indicates that the vector is represented in geographic (or vehicle) frame. The vector
24 [vn , vd ]⊤ is the Earth-relative velocity of point P , represented in the instantaneous tangent plane.
The second term on the right side of (76) is due to the rotation rate of the geographic frame
26 with respect to the Earth-centered inertial frame (that is, transport rate). Equations (75)-(77) are
an example of the nonlinear kinematic model of (1) in the main body of this article.

28 For INS computations, the acceleration vector and angular rate are computed from the

41
 IMU measurements. The IMU consists of a dual-axis accelerometer and a single-axis gyro with
2 output measurements modeled as
⃗u˜v1 = (⃗aviv − ⃗g v ) − ⃗εv1 = f⃗v − ⃗εv1 (78)
ũ2 = ωiv − ε2 , (79)
where the tilde indicates a measurement, f⃗ = ⃗aiv −⃗g is the specific force vector, and ⃗g represents
the gravity vector. The terms
⃗εv1 = −dv1 (⃗u) − z1v and ε2 = −d2 (⃗u) − z2
represent the sum of the deterministic and stochastic accelerometer and gyro errors, as defined
4 in (3).

The IMU provides the measurements in v-frame, but they are needed for computations in
g-frame. Vectors are transformed between frames using a direction cosine matrix (for example,
f⃗g = Rvg f⃗v , where f⃗g is the specific force vector represented in geographic frame). The rotation
matrix Rvg from platform to geographic frame is defined as
" #
cos(θ) sin(θ)
Rvg =
− sin(θ) cos(θ)
and Rgv = (Rvg )⊤ .

6 A navigation system calculates the vehicle state by integration of

" ˙ # " 1 #" #
ϕ̂ 0 v̂n
˙ = R+ĥ (80)
ĥ 0 −1 v̂d
" # " v̂ v̂ #
˙v̂n   n d
g ˜v ˆv ˆg R+ĥ
= R̂v ⃗u1 + ⃗ε1 + ⃗g (ĥ) + −v̂2 (81)
v̂˙ d n
R+ĥ

˙ v̂n
θ̂ = + (ũ2 + ε̂2 ), (82)
R + ĥ
where ε̂1 is an estimate of the accelerometer error vector and ε̂2 is an estimate of the gyro
8 error. Equations (80)-(82) provide an example of (2) in the main body of the article. The INS
integrates these nonlinear equations to propagate the vehicle state through time.

The inputs to (1) and (2) are, respectively,

⊤ h i⊤
⃗u = (⃗uv1 )⊤ , u2 ⃗uˆ = (⃗uˆv1 )⊤ , û2 .

and
Because the actual inputs ⃗uv1 = ⃗aviv − ⃗g v and u2 = ωiv are not available, for the purpose of
integrating (80)-(82), they are computed from the measurements as
⃗uˆv1 = ⃗u˜v1 + ⃗εˆv1 and û2 = ũ2 + ε̂2 .

42
 For use later, define
 
⃗ˆ
fn  = R̂vg ⃗uˆv1 . (83)
ˆ
f⃗d

The calibration terms ⃗εˆv1 and ε̂2 are computed using IMU error model parameters that are esti-
2 mated in real-time. These calibration parameters are denoted as ⃗xd (t) for the IMU deterministic
errors and ⃗xz (t) for the IMU stochastic errors in the augmented state vector defined in (5).

In this example, the vehicle state and estimated vehicle state are defined as
⃗xv = [ϕ, h, vn , vd , θ]⊤ and ⃗xˆv = [ϕ̂, ĥ, v̂n , v̂d , θ̂]⊤ .
Because the actual state ⃗xv is not known, the navigation error state defined as δ⃗x = ⃗xv − ⃗xˆv is
also not known. Defining δ⃗u = ⃗u − ⃗uˆ, the linearized error model of (4) would have
 
0 0 0
0 0 0 
 

 
G(t) =   cos(θ̂) sin( θ̂) 0  , and

− sin( θ̂) cos( θ̂) 0
 
 
0 0 1
 
v̂n 1
 0 − 0 0 
 (R+ĥ)2 R+ĥ 
 

 0 0 0 −1 0 
 ˆ 
 0 − v̂n v̂d v̂d v̂n
f⃗d 
F (t) =  .
 
(R+ĥ)2 ! R+ĥ R+ĥ
 
 2
v̂n ∂⃗g g (h) −2v̂n ˆ 
⃗
 0 + ∂h
0 −f n 
h=ĥ
 
 (R+ĥ)2 R+ĥ 
 
v̂n 1
0 − 0 0
 
(R+ĥ)2 R+ĥ

4 The main purpose of this tutorial is to discuss the issues and methods related to defining the
stochastic error state vector ⃗xz (t) and its state-space model in the form of (6). When this is
6 done, the F and G matrices defined above are used with state augmentation methods to define
the complete error model.

8 The augmented state-space model communicates to the mathematics of the state estimation
process both how the IMU calibration states ⃗xd and ⃗xz change dynamically with time and how
10 they affect the vehicle state estimate ⃗xv . The explicit method by which this is done is by using
(86) to propagate the error state covariance matrix through time. It is critical to note that the
12 vehicle state vector is propagated through time using the nonlinear kinematic model of (2); the
linearized model of (86) is only used to propagate the error covariance matrix.

43
 For more extensions to this example (such as its use to decouple the horizontal and
2 vertical error dynamics to explain the vertical channel instability and horizontal channel Schuler
oscillation), see p. 105 in [3].

4 In actual three-dimensional applications, the attitude representation and its update become
more complicated; nevertheless, the basic approach and issues remain the same.

44
 Sidebar: State Augmentation

2 The estimation algorithm estimates the augmented error state as defined in (5). The
dimension of this augmented state is nx = nv + nd + nz .

4 The state-space error model for the vehicle error state is defined in (4). The state-space
model for the inertial measurement unit (IMU) stochastic errors is defined in (6)-(7). For this
6 sidebar (to allow for modeling three accelerometers and three gyros), the dimension of the output
matrix will change to Cz ∈ ℜ6×nz .

8 Similarly, define the state-space model for the IMU deterministic errors as

⃗x˙ d (t) = Ad ⃗xd (t) + Bd ω

⃗ d (t), (84)
zd (t) = Cd ⃗xd (t), (85)

where Ad ∈ ℜnd ×nd , Bd ∈ ℜnd ×r , and Cd ∈ ℜ6×nd . The parameter r represents the number of
10 distinct noise processes in the deterministic error model. The parameter nd represents the number
of states selected to model the IMU deterministic errors. The elements of ⃗xd in the deterministic
12 error model are usually considered to be unknown constants; therefore, the corresponding model
has Ad , Bd , and r all being identically zero [that is, ⃗x˙ d (t) = 0].

14 Combining (4), (6)-(7), and (84)-(85) the linearized state-space error model is
    
F G Cd G Cz δ⃗xv G 0 " #
⃗
η z
⃗x˙ =  0 0 0   ⃗xd  +  0 0  , (86)
    
ω
⃗z
0 0 Az ⃗xz 0 Bz
where the time dependence of all quantities has been dropped from the notation.

16 Given a set of aiding measurements, the objective of the data fusion system is to estimate
the augmented error state vector x(t) that is defined in (5) in real-time. Success requires that
18 the state vector be observable, which is a well-studied problem [76], [77], [78].

45
 Sidebar: Power Spectral Density

For a stationary process, the correlation function R(τ ) = E⟨x(t)x(t + τ )⟩ and two-sided
PSD are Fourier transform pairs related by
Z ∞
S(ω) = e−jωτ R(τ )dτ (87)
−∞
Z ∞
1
R(τ ) = ejωτ S(ω)dω. (88)
2π −∞
2 Instrument error models include nonstationary stochastic processes such as random walk and
integrated random walk. Nonstationary stochastic processes can be analyzed using average
4 correlation functions and average power spectrum (see Section 2.7 in [11] or p. 109 in [86]),
which are related to each other in the same way as shown in (87)-(88). This article will not
6 distinguish between the two.

46
 Sidebar: Finite-Dimensional Linear State-Space Systems have Even Power Spectra

2 The main text stated that the power spectrum for a linear state-space model (without pure
delay) will be an even polynomial function of s = jω. This sidebar discusses two aspects of
4 this statement.

State-Space to Transfer Function

6 Consider the single-input, single-output, finite-dimensional, linear state-space model

⃗x˙ (t) = F ⃗x(t) + G u(t), (89)
z(t) = H ⃗x(t), (90)
where F ∈ ℜn×n , G ∈ ℜn×1 , and H ∈ ℜ1×n . The parameter n represents the order of the
8 system. The transfer function from u(t) to z(t) is denoted by UZ(s)
(s)
= T (s), where s is the
Laplace variable, Z(s) and U (s) are the Laplace Transforms of z(t) and u(t), respectively. The
10 transfer function T (s) can be computed from the state-space model parameters as (see Section
3.5.2: in [3]):
T (s) = H (sI − F )−1 G. (91)
This transfer function is the ratio of polynomials in s, namely:
N (s)
T (s) = .
D(s)
The purpose of this sidebar is to provide examples to demonstrate that the power spectrum
N (s) N (s∗ )
S(ω) = T (s)T (s∗ )|s=jω =
D(s) D(s∗ ) s=jω
12 is an even polynomial function of ω. For related theory and additional examples, see Sections
3.2-3.7 in [56].

14 State-Space to Power Spectrum

Consider the double integrator state-space system:

" # " #" # " #
ṗ(t) 0 1 p(t) 0
= + u(t)
v̇(t) 0 0 v(t) 1
" #
h i p(t)
z(t) = 1 0 .
v(t)
16 Using eqn. (91), the transfer function is
i s −1 −1 0
" # " #
Z(s) h 1
= 1 0 = 2.
U (s) 0 s 1 s

47
 1
For T (s) = s2
,
1 1 1
S(ω) = = ,
(jω)2 (−jω)2 ω4
which is an even polynomial function of ω. Additional examples can be found in many text
2 books.

Power Spectrum to State-Space

4 The fact that any power spectrum S(ω) that is an even (finite-order) polynomial function of
ω can be represented by a finite-dimensional linear state-space system is shown by first factoring
6 S(ω) = T (jω)T (−jω), where T is the ratio of finite-dimensional polynomials in jω and then
finding a state-space representation for T (s).
2
8 For example, the power spectrum S(ω) = Aω4 can be factored as S(ω) = (jω) A A
2 (−jω)2 , where

the first term provides the transfer T (s) = sA2 . The transfer function UZ(s)
(s)
= sA2 is equivalent to
10 s2 Z(s) = A U (s). Multiplication by s in the Laplace domain corresponds to differentiation in
the time domain. Therefore, z̈(t) = A u(t), which has the state-space model:
" # " #" # " #
ṗ(t) 0 1 p(t) 0
= + u(t)
v̇(t) 0 0 v(t) A
" #
h i p(t)
z(t) = 1 0 .
v(t)

48
 Sidebar: A Brief Historical Review of the Allan Variance

2 The Allan Variance (AV) was originally proposed in the 1960s for the study of frequency
stability of oscillators and signal generators [49], [50], [87]. Since its definition, the AV has
4 found utility for the specification of inertial measurement unit (IMU) performance. The IEEE
standards are written in terms of the AV [28], [29], [30]. This section provides a very brief
6 introduction to the history of the AV.

Consider a signal generator with instantaneous output voltage V (t) given by

V (t) = [V0 + ϵ(t)] sin [2πν0 t + φ(t)] , (92)

8 where V0 and ν0 are the nominal output amplitude and frequency, respectively; and ϵ(t) and φ(t)
are the instantaneous random amplitude and phase fluctuations. From φ(t), the instantaneous
10 fractional frequency fluctuation u(t) is defined as
φ̇(t)
u(t) = . (93)
2πν0
Before the introduction of the AV, the standard measure of frequency stability was the spectral
12 density Su (f ). The AV is an alternative, time-domain measure of frequency stability defined as
* N N
!2 +
1 X 1 X
σu2 (τ : N, Ts ) = ūk − ūj . (94)
N − 1 k=1 N j=1

The notation σu2 (N, Ts , τ ) is standard in the AV literature. This sidebar uses the more descriptive
14 notation σu2 (τ : N, Ts ) to indicate that N and Ts are parameters that the analyst selects to evaluate
the value at cluster duration τ . In this notation, N is the number of clusters of duration τ that
16 are used in the computation, and Ts is the time between the start of consecutive clusters. The
operator ⟨·⟩ indicates an infinite time-average. In (94),
tZ
k +τ
1
ūk = u(t)dt, (95)
τ
tk

18 with tk = tk−1 + Ts . The expression in (94) can be understood as the sample variance of N
averages of u(t) each over a time interval of duration τ . The minimum value of τ and Ts is the
20 sample period T . The symbols used in this discussion are defined in Table 3. The relationship
between these parameters is illustrated in Figure 12. In the case where the time interval satisfies
22 Ts > τ , the computation has deadtime between clusters where data is unused.

The “Two-Sample Without Dead-Time” formula

* +
2
(ūk+1 − ūk )
σu2 (τ ) ≡ σu2 (τ : 2, τ ) = (96)
2

49
 (that is, ⟨σu2 (N = 2, Ts = τ, τ )⟩ in the standard notation) is recommended in [49], [50] and
2 eventually became known as the AV.

Since practical data records are of finite length, the infinite time averages are not available;
4 therefore, approximations are required. Barnes et al. [50] proposed and studied
m−1
1 X
σ̂u2 (τ ) = (ūk+1 − ūk )2 . (97)
2(m − 1) k=1
Note the hat in the left side of (97) which indicates that it is designed as an estimate of σu2 (τ )
6 defined in (96). This formula has become known as the “Non-Overlapping” AV, in reference to
the fact that Ts = τ .

8 Because the duration of the available data set is finite, for each value of τ and choice of
Ts , the number of clusters m will change (with longer clusters and larger values of Ts yielding
10 smaller values of m). Given a data record with L samples and a constant sample period T , the
total experiment duration is L T seconds. The number of averages ūk (that is, clusters) that can
12 be computed for a cluster duration of τ = nT , without deadtime, is m = L/n.

Several alternative AV formulae have been proposed based on different choices for Ts . For
14 instance, Howe, Allan, and Barnes [88] introduced the (Fully) Overlapping AV with Ts = T
[which is the same as (12)]:
L−2n
1 X
σ̂u2 (τ ) = (ūk+n − ūk )2 . (98)
2(L − 2n) k=1
16 Its stated objective is to provide the best confidence in the estimates, which is achieved by high
data utilization, that is, the number of formed averages ūk is no longer m = L/n, but instead,
18 m = L − n; therefore, the estimation accuracy of the Overlapping AV (based on (L − 2n)
differences) increases dramatically (relative to the Non-Overlapping AV) for long cluster-times.
20 The Overlapping AV has since become the standard for IMU stochastic error modeling [28],
[29].

22 Allan and Barnes [87] modified the Overlapping AV to improve its ability to distinguish
stochastic processes with spectral densities Su (f ) ≈ f α , such that α = +1 or +2 (flicker phase
24 noise and white phase noise, respectively). However, these are generally not the main sources
of stochastic errors corrupting an IMU [see (17)].

26 A more recent alternative to (98) is the Not-Fully-Overlapping AV, wherein τ > Ts > T
[54]. In the article, the authors show that the method has similar estimation accuracy to the
28 (Fully) Overlapping AV (but at a reduced computational cost), which is relevant because AV
analysis for IMU characterization generally requires large datasets.

50
 It is not the purpose of this tutorial to detail the historical development of the AV [48],
2 [89], [90]. For the interested reader, Table 4 summarizes and compares various formulations of
the AV and provides pointers to selected references.

51
 Sidebar: Discussion of Equation (62)

2 Equation (62) may appear to be counterintuitive. Why would the discrete-time measurement
variance Qηd decrease as the sample period T increases? This phenomenon has a long history
4 that can be understood from different perspectives.

Sensor Model. The general assumption is that a discrete-time measurement is obtained as the
mean of the continuous-time quantity within the sample interval:
1 tk+1
Z
ũ(k) = (u(τ ) + ηω (τ )) dτ (99)
T tk
1 tk+1
Z
ũ(k) = ū(k) + ηω (τ )dτ (100)
T tk
ũ(k) = ū(k) + η(k), (101)

where the discrete-time measurement noise is

1 tk+1
Z
η(k) = ηω (τ )dτ. (102)
T tk
If ηω (τ ) is white, then its covariance function is E⟨ηω (ζ)ηω (τ )⟩ = SN δ(ζ − τ ), where δ denotes
the Dirac delta function and SN is the power spectral density (PSD). Therefore, the covariance
of η(k) is computed as

Qηd (k) = E ⟨η(k)η(k)⟩

 Z tk+1   Z tk+1 
1 1
=E ηω (τ )dτ ηω (ζ)dζ
T tk T tk
Z tk+1 Z tk+1
1
= 2 E ⟨ηω (τ )ηω (ζ)⟩ dζdτ
T tk tk
Z tk+1 Z tk+1
1
= 2 SN δ(ζ − τ )dζdτ
T tk tk
Z tk+1
1
= 2 SN dτ
T tk
1
Qηd (k) = SN , (103)
T
which is the same as (62).

6 Angle Increments. From (99), the discrete-time samples from the inertial measurement unit
(IMU) may be presented (that is, scaled) as either an angular rate (or acceleration) measurement
8 u(k) or an angle (or velocity) increment ∆(k) = u(k) T over the time increment of length T .
The analysis of (56)-(62) presented the IMU white noise conversion from continuous to discrete

52
 time for the first case: IMU angular rate (or acceleration) outputs, which resulted in (62). This
2 section considers the analysis for the case where IMU outputs angle (or velocity) increments.

The discrete-time model that is equivalent to (51) is

xv (k + 1) = xv (k) + ∆(k). (104)

4 The discrete-time model that is equivalent to (52) is

˜
x̂v (k + 1) = x̂v (k) + ∆(k). (105)

Scaling both sides of (58) by T yields the discrete-time measurement model,

˜
∆(k) = ∆(k) + η∆ (k), (106)

6 with white measurement noise η∆ (k) ∼ N (0, Qη∆ ). By the definition of ∆(k) in the previous
paragraph, Qη∆ = T 2 Qηd , where Qηd is defined in (45). The error signal e(k) = xv (k) − x̂v (k)
8 has the time propagation model

e(k + 1) = e(k) − η∆ (k) (107)

which is a discrete-time random walk process. The discrete-time propagation of the covariance
10 of e(k) driven by η∆ (k) is

Pe (k + 1) = Pe (k) + Qη∆ for any k ≥ 0. (108)

Due to the assumption that the initial covariance of e(k) is zero, (108) simplifes to

Pe (k) = k Qη∆ . (109)

12 Because the continuous and discrete-time models are equivalent, their covariance must be
the same at the discrete sample times. Equating (55) to (109) yields

Pe (t) t=kT
= Pe (k)
SN k T = k Qη∆ ,

14 which provides
Qη∆ = SN T, (110)

which is equivalent to (62) because Qη∆ = T 2 Qηd .

Q
16 The fact that the PSD SN must be equal to Tη∆ is discussed in Example 3.20 in [91],
which attributes that example to Kalman in [92]. The example discusses continuous-time white
18 noise as the limit of discrete-time white noise as T approaches zero.

Unit Analysis. Consider the units of SN , Qη∆ , and Qηd .

53
 2 2
• The symbol SN represents the PSD of ηz (t), which has units of (deg/s)
Hz
= (deg)
s
for gyros
2
(m/s )
2 2
2 and Hz = (m) s3
for accelerometers.
• The symbol Qηd represents the covariance of η(k), which has units of (deg/s)2 for gyros
2
4 and (m/s2 ) for accelerometers.
• The symbol Qη∆ represents the covariance of η∆ (k), which has units of (deg)2 for gyros
6 and (m/s)2 for accelerometers.

Note that all these units work out consistently in (62) and (110). Equation (62) is used to compute
8 the covariance of the discrete-time white noise covariance Qηd (which is needed for the design
of the state estimator) from the continuous-time PSD SN [which is extracted from the Allan
10 Standard Deviation (ASD)].

54
 Short Biography for each Author

2 Farrell

Jay A. Farrell received B.S. degrees in physics and electrical engineering from Iowa
4 State University, and M.S. and Ph.D. degrees in electrical engineering from the University
of Notre Dame. At Charles Stark Draper Lab (1989-1994), he was principal investigator
6 on projects involving autonomous vehicles, receiving the Engineering Vice President’s Best
Technical Publication Award in 1990 and Recognition Awards for Outstanding Performance and
8 Achievement in 1991 and 1993. He is the KA Endowed Professor in the Department of Electrical
and Computer Engineering at the University of California, Riverside. For the IEEE Control
10 Systems Society (CSS), he has served as vice president finance, vice president of technical
activities, CSS general vice chair of IEEE CDC-ECC 2011, general chair of IEEE CDC 2012,
12 president elect, president, and past president. For IEEE he served three terms on the Fellow
Committee, as EAB treasurer, and on IEEE FinComm. In 2020-2021, he served as president of
14 AACC. He is author of over 250 technical articles and three books. He was recognized as a
GNSS Leader to Watch by GPS World Magazine in 2009 and is a Distinguished Member of
16 IEEE CSS, a Fellow of the IEEE, a Fellow of AAAS, and a Fellow of IFAC.

Silva

18 Felipe O. Silva received the B.S. Degree in automatic control engineering (with honors)
from the Federal University of Itajubá, in Brazil (2012); the M.S. degree in systems engineering
20 from the National Institute of Applied Sciences, Center Val de Loire, in France (2013); and
the Ph.D. degree in aeronautical and mechanical engineering from the Aeronautics Institute of
22 Technology, Brazil (2016). From 2013 to 2014, he was an assistant researcher with the Institute
of Aeronautics and Space, in Brazil. Since 2014, he has been an assistant professor with the
24 Federal University of Lavras, in Brazil. Since 2017, he has been an associate with the technology-
based microenterprise Inovação em Mecanização Agrı́cola Ceifa Ltda., in Brazil. In 2018, he
26 was a visiting professor with the Central School of Nantes, in France. From 2018 to 2019, he
was a visiting professor with the University of California Riverside, in Riverside, USA. His
28 research area includes: guidance, navigation and control systems; inertial navigation systems;
global navigation satellite systems; sensor fusion; instrumentation; robotics; hydro-pneumatic
30 systems and precision agriculture. Since 2020, he has been a member of the Brazilian Computer
Society.

55
 Rahman

2 Farzana Rahman received the Ph.D. degree in electrical and computer engineering from
University of California, Riverside, in 2020 and the B.S. degree in electrical engineering from
4 Bangladesh University of Engineering and Technology in 2014. She currently works as a software
engineer at Nuro, working on robotics and self-driving vehicles. Her current research interests
6 include robot localization and mapping algorithm design.

Wendel

8 Jan Wendel received the Dipl.-Ing. and Dr.-Ing. degrees in electrical engineering from the
University of Karlsruhe, Germany in 1998 and 2003, respectively. From 2003 to 2006, he was
10 an assistant professor at the University of Karlsruhe, where he was private lecturer until 2018.
In 2006, he joined MBDA in Munich, Germany, working on the development of navigation
12 algorithms for various applications. In 2009, he joined Airbus Defense and Space in Munich,
where he is involved in various activities related to satellite navigation, including acquisition and
14 tracking algorithms, interference detection and characterization, integrity algorithms, and PRS
receiver development.

56
 Article Summary

2 Autonomous vehicle technology is advancing rapidly. Their control capabilities often

rely on high-bandwidth state estimation incorporating inertial measurements. High-performance
4 state estimation incorporates inertial measurement error models through the process of state
augmentation to enable on-the-fly instrument calibration.

6 This article is a tutorial describing the process and issues related to developing a state-
space model for the stochastic errors affecting an Inertial Measurement Unit (IMU). The starting
8 point is the instrument error characterization data sheet provided by the manufacturer, which is
typically either a Allan Standard Deviation graph or the Allan Variance parameters extracted
10 from that graph. The desired output of the modeling process is a linear discrete-time state-space
model of the IMU stochastic errors suitable for augmentation to the INS error state model.

12 Along with this tutorial, supplementary open source software is available. One software
component does the following: (1) Given a continuous-time state-space IMU stochastic error
14 model selected by the designer to match the Allan variance, the software computes a discrete-
time equivalent state-space model. (2) Given that discrete-time model, it produces a stochastic
16 error sequence suitable for Allan Variance computations. (3) Given a sequence of stochastic
errors, it computes and plots the Allan Variance. Given Allan Variance data and a specific
18 continuous-time state-space IMU stochastic error model structure, a second software component
implements an optimization-based approach to select the model parameters to match the Allan
20 Variance data.

57
TABLE 1: Dominant errors in consumer grade inertial measurement units (IMUs)

Allan Variance PSD

Noise type Coef.
Acc: (m/s2 )2 Acc: m2 /s3
(Coef.) unit
Gyro: (deg/s)2 Gyro: deg2 /s
Ang./Vel. Acc: m/s3/2 N2
N2
random walk, N Gyro: deg/s1/2 τ
Bias Acc: m/s2 2 B 2 ln(2) B2
instability, B Gyro: deg/s π 2πf
2
Acc: m/s5/2

Rate/Accel. K2 τ K
random walk, K Gyro: deg/s3/2 3 2πf

58
TABLE 2: Extracted Allan Variance (AV) parameters related to Figure 8

Coef. N B K TB
Untuned Value 0.0033 0.0011 0.00014 32
Manually-tuned Value 0.0033 0.0004 0.00014 20
Optimization-based Value 0.0033 0.0001 0.00012 50
Unit m/s3/2 m/s2 m/s5/2 s

59
TABLE 3: Symbols used in the discussion of Allan Variance (AV).

Symbol Meaning
T Sample period for u(t)
L Total number of samples in the data set
LT Duration of the data set
τ Cluster or averaging duration
n Number of sample periods per cluster: τ = nT
N Number of clusters of duration τ used in (94)
m Number of values of σu2 (τ : 2, τ ) averaged in (97)
tk Start time of the k-th averaging interval, see (95)
Ts Time interval between consecutive averaging intervals
ūk Average of u(t) for t ∈ [tk , tk + τ ], see (95)

60
 TABLE 4: Main Allan Variance (AV) estimation formulae

AV Stated Selected
Estimate Benefits References
Non-Overlapping
[49], [38],
with Dead- —
[93], [94], [95], [96], [97]
Time (NODT)
Non-Overlapping Simpler to compute [50], [51], [52]
(NO) (compared to NODT) [98], [94], [99]
Gives high
data utilization/
Overlapping
better confidence [88], [100], [101], [102]
(O)
in the estimate
(compared to NO)
Able to better
Modified distinguish some [87], [102]
(M) types of noise [103]
(compared to O)
Not-Fully- More computational
Overlapping efficient [54]
(NFO) (compared to O)

61
Figure 1: Allan Standard Deviation (ASD) plots for three Crossbow µNav gyros [36], [45], with
straight line approximations for dominant errors. Sampling interval is T = 0.02 s.

62
Figure 2: Accelerometer Allan Standard Deviation (ASD) plot for the NV IMU-1000 from Nav
Technology, with straight line approximations for dominant errors. Sampling interval is T = 0.01
s.

63
Figure 3: Typical gyro Allan Standard Deviation (ASD) shape corresponding to (17)-(19).

64
Figure 4: Angular/Velocity random walk Allan Standard Deviation (ASD) plot. See (23).

65
Figure 5: Rate/Acceleration random walk Allan Standard Deviation (ASD) plot. See (27).

66
Figure 6: Bias instability Allan Standard Deviation (ASD) plot with f0 = 1. See (31).

67
 √
Figure 7: Gauss Markov Allan Standard Deviation (ASD) plot with qB = SB TB . See (35).

68
Figure 8: Allan Standard Deviation (ASD) from Figure 2 along with ASD plots for N, B, and
K.

69
Figure 9: Allan Standard Deviation (ASD) plot for inertial measurement unit (IMU) data, model
defined by (42), and simulated data as described in the “ASD Verification of the State-Space
Model” section.

70
Figure 10: Generic block diagram of a global navigation satellite system (GNSS)-aided inertial
navigation system (INS).

71
 q
u

w P
n

h
d
R
z
f
p

Figure 11: Variables for a two-dimensional simplified inertial navigation system (INS) example.

72
Figure 12: Relationship between T , τ , Ts , tk , and ūk for this example: τ = 5T and Ts = 7T ;
therefore, it has a deadtime of 2T .

73