# Customizable Stochastic High Fidelity Model of the Sensors and Camera
onboard a Low SWaP Fixed Wing Autonomous Aircraft

Eduardo Gallo111The author holds a MSc in Aerospace Engineering by the
Polytechnic University of Madrid and has twenty-two years of experience
working in aircraft trajectory prediction, modeling, and flight simulation. He
is currently a Senior Trajectory Prediction and Aircraft Performance Engineer
at Boeing Research & Technology Europe (BR&TE), although he is publishing this
article in his individual capacity and time as part of his PhD thesis titled
“Autonomous Unmanned Air Vehicle GNSS-Denied Navigation”, advised by Dr.
Antonio Barrientos within the Centre for Automation and Robotics of the
Polytechnic University of Madrid. 222Contact<EMAIL_ADDRESS>https://orcid.org/0000-0002-7397-0425

(February 2021)

## Abstract

The navigation systems of autonomous aircraft rely on the readings provided by
a suite of onboard sensors to estimate the aircraft state. In the case of
fixed wing vehicles, the sensor suite is composed by triads of accelerometers,
gyroscopes, and magnetometers, a Global Navigation Satellite System (GNSS)
receiver, and an air data system (Pitot tube, air vanes, thermometer, and
barometer), and is often complemented by one or more digital cameras. An
accurate representation of the behavior and error sources of each of these
sensors, together with the images generated by the cameras, in indispensable
for flight simulation and the evaluation of novel inertial or visual
navigation algorithms, and more so in the case of low SWaP (size, weight, and
power) aircraft, in which the quality and price of the sensors is limited.
This article presents realistic and customizable models for each of these
sensors, which have been implemented as an open-source
$\displaystyle\mathrm{{C\nolinebreak[4]\raisebox{1.72218pt}{\tiny\bf++}}}$
simulation, available at [1]. Provided with the true variation of the aircraft
state with time, the simulation provides a time stamped series of the errors
generated by all sensors, as well as realistic images of the Earth surface
that resemble those taken from a real camera flying along the indicated state
positions and attitudes.

##  1 Introduction and Outline

The sensors onboard an autonomous aircraft measure various aspects of the
aircraft _real_ or _actual trajectory_
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
TRUTH}}$ and provide these measurements to the aircraft guidance, navigation
and control (GNC) system. The outputs of these sensors, collectively known as
the _sensed trajectory_
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}}=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
SENSED}}$, represent the only link between the real but unknown actual
trajectory and the GNC system in charge of achieving an actual trajectory that
deviates as little as possible from the guidance targets.

SENSORS$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}\left(t_{t}\right)=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
TRUTH}\left(t_{t}\right)}$$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}}\left(t_{s}\right)=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
SENSED}\left(t_{s}\right)}$ Figure 1: Sensors flow diagram

Note that although the actual aircraft state varies continuously in the real
world, in simulation it is usually the outcome of a high frequency discrete
integration process [2] that results in
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}\left(t_{t}\right)=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
TRUTH}\left(t_{t}\right)}$, where $\displaystyle\mathrm{t_{t}=t\cdot\Delta
t_{\scriptscriptstyle TRUTH}}$. The sensed trajectory
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}}\left(t_{s}\right)=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
SENSED}\left(t_{s}\right)}$, where $\displaystyle\mathrm{t_{s}=s\cdot\Delta
t_{\scriptscriptstyle SENSED}}$, is however intrinsically discrete, although
the working frequency of the different sensors may vary. This article
considers that all sensors operate at the same rate of
$\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$, with the
exception of the GNSS receiver and the onboard camera, which work at
$\displaystyle\mathrm{\Delta t_{\scriptscriptstyle GNSS}}$ and
$\displaystyle\mathrm{\Delta t_{\scriptscriptstyle IMG}}$, respectively. It is
also assumed that all sensors are fixed to the aircraft structure in a
strapdown configuration and that their measurement processes are instantaneous
and time synchronized with each other at their respective frequencies. Table 1
provides the default frequencies employed in this article as well as in the
corresponding software implementation [1], where all values can be adjusted by
the user.

Discrete Time | Frequency | Rate
---|---|---
$\displaystyle\mathrm{t_{t}=t\cdot\Delta t_{\scriptscriptstyle TRUTH}}$ | $\displaystyle\mathrm{500\left[hz\right]}$ | $\displaystyle\mathrm{0.002\left[sec\right]}$
$\displaystyle\mathrm{t_{s}=s\cdot\Delta t_{\scriptscriptstyle SENSED}}$ | $\displaystyle\mathrm{100\left[hz\right]}$ | $\displaystyle\mathrm{0.01\left[sec\right]}$
$\displaystyle\mathrm{t_{i}=i\cdot\Delta t_{\scriptscriptstyle IMG}}$ | $\displaystyle\mathrm{10\left[hz\right]}$ | $\displaystyle\mathrm{0.1\left[sec\right]}$
$\displaystyle\mathrm{t_{g}=g\cdot\Delta t_{\scriptscriptstyle GNSS}}$ | $\displaystyle\mathrm{1\left[hz\right]}$ | $\displaystyle\mathrm{1\left[sec\right]}$
Table 1: Default frequencies of the states in the actual and sensed
trajectories

The sensed trajectory can be defined as a time stamped series of state vectors
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}}=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
SENSED}}$ that groups the measurements provided by the different onboard
sensors, comprising the only view of the actual trajectory at the disposal of
the GNC system:

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}}=\mbox{\boldmath{$\displaystyle\mathrm{x}$}}_{\scriptscriptstyle
SENSED}=\left[\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B},\,\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B},\,\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}}^{\scriptscriptstyle
B},\,\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}}_{\scriptscriptstyle
GDT},\,\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}}^{\scriptscriptstyle
N},\,\widetilde{p},\,\widetilde{T},\,\widetilde{v}_{\scriptscriptstyle
TAS},\,\widetilde{\alpha},\,\widetilde{\beta},\,\,\mbox{\boldmath{$\displaystyle\mathrm{I}$}}\right]^{T}}$
(1)

Note that when present, the super index represents the frame or reference
system in which a certain variable is viewed; if two sub indexes are present,
it implies that the vector goes from the first frame to the second. Table 1
makes use of the body frame $\displaystyle\mathrm{F_{\scriptscriptstyle B}}$,
rigidly attached to the aircraft structure with origin in its center of mass
[3], the NED frame $\displaystyle\mathrm{F_{\scriptscriptstyle N}}$ also
centered on the aircraft center of mass with axes in the North-East-Down
directions [3], and the inertial frame
$\displaystyle\mathrm{F_{\scriptscriptstyle I}}$, usually considered as
centered in the Sun with axes fixed with respect to other stars [4].
Additionally, the _specific force_
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle{IB}}}$
is defined as the non gravitational acceleration experienced by the aircraft
body with respect to an inertial frame [5].

Components | Variable | Measured by | Acronym | Rate
---|---|---|---|---
Specific Force | $\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle IB}^{\scriptscriptstyle B}}$ | Accelerometers | ACC | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Inertial Angular Velocity | $\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle IB}^{\scriptscriptstyle B}}$ | Gyroscopes | GYR | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Magnetic Field | $\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}}^{\scriptscriptstyle B}}$ | Magnetometers | MAG | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Geodetic Coordinates | $\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{x}$}}}_{\scriptscriptstyle GDT}}$ | GNSS Receiver | GNSS | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle GNSS}}$
Ground Velocity | $\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}}^{\scriptscriptstyle N}}$ | GNSS Receiver | GNSS | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle GNSS}}$
Air Pressure | $\displaystyle\mathrm{\widetilde{p}}$ | Barometer | OSP333OSP stands for Outside Static Pressure. | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Air Temperature | $\displaystyle\mathrm{\widetilde{T}}$ | Thermometer | OAT444OAT stands for Outside Static Temperature. | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Airspeed | $\displaystyle\mathrm{\widetilde{v}_{\scriptscriptstyle TAS}}$ | Pitot Tube | TAS555TAS stands for True Air Speed. | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Angle of Attack | $\displaystyle\mathrm{\widetilde{\alpha}}$ | Air Vanes | AOA666AOA stands for Angle Of Attack. | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Angle of Sideslip | $\displaystyle\mathrm{\widetilde{\beta}}$ | Air Vanes | AOS777AOS stands for Angle of Sideslip. | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle SENSED}}$
Image | $\displaystyle\mathrm{I}$ | Digital Camera | CAM888CAM stands for Camera. | $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle IMG}}$
Table 2: Components of sensed trajectory

The main objective of this article is to properly model the sources of error
incurred by the different sensors onboard the aircraft, with special emphasis
on the inertial sensors (accelerometers and gyroscopes) that provide the
essential information on which the navigation system of low SWaP vehicles
relies. The model parameters need to be _customizable_ so the user can employ
the values that better resemble the performances of the specific equipment
being modeled. They also need to be _stochastic_ to properly represent the
nature of the different random processes involved, while ensuring that the
time variation of the errors generated by each sensor can be repeated if so
desired. Finally, the primary use of the sensor models developed in this
article is to be part of a _high fidelity_ flight simulation such as that
described in [2], so the models need to be realistic and include as few
simplifications as possible.

Section 2 describes the objectives of the simulation with special emphasis on
its stochastic properties and how it can be customized for different sensors.
It is followed by detailed explanations of the models representing the errors
present in the measurements of the different sensors. Section 3 focuses on the
inertial sensors and explains their characteristics and error sources,
resulting in a model that captures their stochastic properties. A similar but
less detailed approach is taken in section 4 with the magnetometers, air data
system, and GNSS receiver. Section 5 describes the characteristics of the
digital camera employed for visual navigation, and the tool employed to
generate realistic images that resemble what a real camera would view if
located at the same position and attitude; although the camera differs from
all other sensors in that it does not provide a measurement or reading but a
digital image, in this article it is indeed considered a sensor as it provides
the GNC system with information about its surroundings that can be employed
for navigation. Last, section 6 describes sensor calibration activities that
are indispensable for the determination of various of the parameters present
in the sensor models.

##  2 Customizable and Stochastic Simulation

The different sensor models described in this article have been implemented in
a C++ simulation (the open-source code is available at [1]) with the objective
of properly simulating their stochastic nature while ensuring that the
simulated variation with time of the errors can be repeated at a later time if
so desired. The implementation hence follows these steps:

  1. 1.

Initialize a discrete uniform distribution with any seed (any value is valid,
so 1 was employed by the author), which produces random integers where each
possible value has an equal likelihood of being produced. Call this
distribution a number of times equal or higher than twice the maximum number
of runs to be executed (each run provides the variation of time of all sensors
for an unlimited amount of time, and corresponds to a single aircraft flight),
divide them into two groups of the same size, and store the results for later
use. These values, called $\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle
A}}$ and $\displaystyle\mathrm{\Upsilon_{j,\scriptscriptstyle F}}$, are
respectively the _aircraft seeds_ and the _flight seeds_ , where i is the
aircraft number representing given fixed error realizations999Fixed error
contributions vary from aircraft to aircraft but are constant for all flights
of that aircraft., and j is the trajectory number representing run-to-run and
in-run error realizations101010Run-to-run and in-run error contributions vary
from one flight of a given aircraft to the next.. The stored aircraft and
flight seeds become the initialization seeds for each of the simulation
executions or runs, so this step does not need to be repeated.

  2. 2.

Every time the simulator needs to obtain the errors generated by the different
sensors (which usually correspond to a given flight), it is initialized with a
given aircraft seed $\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$
together with a flight seed
$\displaystyle\mathrm{\Upsilon_{j,\scriptscriptstyle F}}$. As these seeds are
the only inputs required for all the stochastic processes within the sensors,
the results of a given run can always be repeated by employing the same seeds.

Type |  | Error Sources | Seeds
---|---|---|---
Aircraft | i | fixed | $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle ACC},\,\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle GYR},\,\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle MAG},\,\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle PLAT},\,\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle CAM}}$
Flight | j | run-to-run & in-run | $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle ACC},\,\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle GYR},\,\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle MAG},\,\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle OSP},\,\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle OAT}}$
$\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle
TAS},\,\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle
AOA},\,\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle
AOS},\,\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle GNSS}}$
Table 3: Sensor seeds

The selected aircraft and flight seeds are then employed to initialize two
different discrete uniform distributions. One is executed five times to
provide the _fixed sensor seeds_
($\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle
ACC}}$, $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle
A,\scriptscriptstyle GYR}}$,
$\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle
MAG}}$, $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle
A,\scriptscriptstyle PLAT}}$, and
$\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle
CAM}}$), while the other is realized nine times to obtain the _run sensor
seeds_ ($\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle
F,\scriptscriptstyle ACC}}$,
$\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle
GYR}}$, $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle
F,\scriptscriptstyle MAG}}$,
$\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle
OSP}}$, $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle
F,\scriptscriptstyle OAT}}$,
$\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle
TAS}}$, $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle
F,\scriptscriptstyle AOA}}$,
$\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle
AOS}}$, and $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle
F,\scriptscriptstyle GNSS}}$). These seeds hence become the initialization
seeds for each of the different sensors described throughout this article.

  3. 3.

Each sensor relies on either one or two standard normal distributions
$\displaystyle\mathrm{N\left(0,\,1\right)}$, depending on whether its error
model is based exclusively on run-to-run and in-run error contributions or it
also contains fixed error sources. The normal distributions of every sensor
are initialized with the corresponding seeds
($\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle F}}$ and
$\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle R}}$) for that sensor.

  4. 4.

Upon initialization, the _fixed_ normal distribution of every sensor is
employed to generate all the values corresponding to scale factors, cross
couplings, hard iron magnetism, and mounting errors. The _run_ normal
distribution in turn is employed to generate the required bias offsets.

  5. 5.

Once the simulation has been initialized, it is able to estimate the errors
generated by each sensor working at the required sensor rate. As the
simulation time advances, every time a sensor is called to provide a
measurement, its already initialized and used run normal distribution is
called to generate the corresponding random walk increments and white noises.

  6. 6.

Additionally, when provided with a time series of the aircraft position and
attitude, the simulation generates realistic images of the Earth surface that
resemble those taken with an onboard camera.

All the error models described in the following sections rely on a series of
parameters whose value varies depending on the specific sensors represented in
the simulation. Even though the user is expected to specify the value of these
parameters, this article suggests the values used by the author to represent
the sensors that could nowadays be installed onboard a small size fixed wing
low SWaP autonomous aircraft, and which from this point on are referred as the
“ _Baseline_ ” configuration.

##  3 Inertial Sensors

Inertial sensors comprise _accelerometers_ and _gyroscopes_ , which measure
specific force and inertial angular velocity about a single axis, respectively
[6]. An _inertial measurement unit_ (IMU) encompasses multiple accelerometers
and gyroscopes, usually three of each, obtaining three dimensional
measurements of the specific force and angular rate [3] viewed in the platform
frame $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ (section 3.5). The
individual accelerometers and gyroscopes however are not aligned with the
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ axes, but with those of the
non orthogonal accelerometers $\displaystyle\mathrm{F_{\scriptscriptstyle A}}$
and gyroscopes $\displaystyle\mathrm{F_{\scriptscriptstyle Y}}$ frames, also
defined in section 3.5. The output of the inertial sensors must hence first be
transformed from the $\displaystyle\mathrm{F_{\scriptscriptstyle A}}$ and
$\displaystyle\mathrm{F_{\scriptscriptstyle Y}}$ frames to
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$, as described in sections 3.6
and 3.7, and then from the $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$
frame to the body frame $\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ as
explained in section 3.8, where they can be employed by the navigation system.
The accelerometers and gyroscopes are assumed to be infinitesimally small and
located at the IMU reference point (section 3.8), which coincides with the
origin of these three frames $\displaystyle\mathrm{\left(O_{\scriptscriptstyle
P}=O_{\scriptscriptstyle A}=O_{\scriptscriptstyle Y}\right)}$.

The IMU is physically attached to the aircraft structure in a strapdown
configuration, so both the displacement
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}$ and the Euler angles
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}^{\scriptscriptstyle
BP}=\left[\psi^{\scriptscriptstyle P},\ \theta^{\scriptscriptstyle P},\
\xi^{\scriptscriptstyle P}\right]^{T}}$ that describe the relative position
and rotation between the body $\displaystyle\mathrm{F_{\scriptscriptstyle B}}$
and platform $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ frames are
constant. Accelerometers can be divided by their underlying technology into
pendulous and vibrating beams, while gyroscopes are classified into spinning
mass, optical (ring laser or fiber optic), and vibratory [5]. Current inertial
sensor development is mostly focused on _micro machined electromechanical
system_ (MEMS) sensors111111There exist both pendulous and vibrating beam MEMS
accelerometers, but all MEMS gyroscopes are vibratory., which makes direct use
of the chemical etching and batch processing techniques used by the
electronics integrated circuit industry to obtain sensors with small size, low
weight, rugged construction, low power consumption, low price, high
reliability, and low maintenance [7]. On the negative side, the accuracy of
MEMS sensors is still low, although tremendous progress has been achieved in
the last two decades and more is expected in the future.

There is no universal classification of inertial sensors according to their
performance, although they can be broadly assigned into five different
categories or grades: marine (submarines and spacecraft), aviation (commercial
and military), intermediate (small aircraft and helicopters), tactical
(unmanned air vehicles and guided weapons), and automotive (consumer) [5]. The
full range of grades covers approximately six orders of magnitude of gyroscope
performance and only three for the accelerometers, but higher performance is
always associated with bigger size, weight, and cost. Tactical grade IMUs
cover a wide range of performance values, but can only provide a stand-alone
navigation solution for a few minutes, while automotive grade IMUs are
unsuitable for navigation.

The different errors that appear in the measurements provided by
accelerometers and gyroscopes are described in section 3.1. Section 3.2
presents a model for the measurements of a single inertial sensor, while
sections 3.3 and 3.4 focus on how to obtain the specific values for white
noise and bias on which the model relies from the documentation. Section 3.5
describes the reference systems required to represent the IMU measurements.
Additional errors appear when three accelerometers or gyroscopes are employed
together, and these are modeled in section 3.6 for accelerometers and section
3.7 for gyroscopes. The analysis of the inertial sensors concludes with
sections 3.9 and 3.10, which provide a comprehensive error model for the IMU
measurements, as well as the specifications of the baseline configuration. The
final model also depends on the relative position of the IMU with respect to
the body frame, which is described in section 3.8.

###  3.1 Inertial Sensor Error Sources

In addition to the accelerometers and gyroscopes, an IMU also contains a
processor, storage for the calibration parameters, one or more temperature
sensors, and a power supply. As described below, each sensor has several error
sources, but each of them has four components: _fixed_ contribution,
_temperature_ dependent variation, _run-to-run_ variation, and _in-run_
variation [5, 8]. The first two can be measured at the laboratory (at
different temperatures) and the calibration results stored in the IMU so the
processor can later compensate the sensor outputs based on the reading
provided by the temperature sensor. Calibration however increases
manufacturing costs so it may be absent in the case of inexpensive sensors.
The run-to-run variation results in a contribution to a given error source
that varies every time the sensor is employed but remains constant within a
given run. It can not be compensated by the IMU processor but can be
calibrated by the navigation system every time it is turned on with a process
known as fine-alignment [3, 5, 9]. The in-run contribution to the error
sources slowly varies during execution and can not be calibrated in the
laboratory nor by the navigation system.

IMU Grade |  Accelerometer Bias $\displaystyle\mathrm{\left[mg\right]}$ |  Gyroscope Bias $\displaystyle\mathrm{\left[deg/hr\right]}$
---|---|---
Marine | 0.01 | 0.001
Aviation | 0.03 - 0.1 | 0.01
Intermediate | 0.1 - 1 | 0.1
Tactical | 1 - 10 | 1 - 100
Automotive | > 10 | > 100
Table 4: Typical inertial sensor biases according to IMU grade

Let’s now discuss the different sources of error that influence an inertial
sensor [5, 10, 11, 12]:

  * •

The _bias_ is an error exhibited by all accelerometers and gyroscopes that is
independent of the underlying specific force or angular rate being measured,
and comprises the dominant contribution to the overall sensor error. It can be
defined as any nonzero output when the sensor input is zero [8], and can be
divided into its static and dynamic components. The static part, also known as
fixed bias, _bias offset_ , turn-on bias, or bias repeatability, comprises the
run-to-run variation, while the dynamic component, known as in-run bias
variation, _bias drift_ , bias instability (or stability), is typically about
ten percent of the static part and slowly varies over periods of order one
minute. As the bias is the main contributor to the overall sensor error, its
value can be understood as a sensor quality measure. Table 3.1 provides
approximate values for the inertial sensor biases according to the IMU grade
[5]:

While the bias offset can be greatly reduced through fine-alignment [3, 5, 9],
the bias drift can not be determined and needs to be modeled as a stochastic
process. It is mostly a warm up effect that should be almost non existent
after a few minutes of operation, and corresponds to the minimum point of the
sensor Allan curve [11, 13]. It is generally modeled as as a random walk
process obtained by the integration of a white noise signal, coupled with
limits that represent the conclusion of the warm up process.

  * •

The _scale factor error_ is the departure of the input output gradient of the
instrument from unity following unit conversion at the IMU processor. It
represents a varying relationship between sensor input and output caused by
aging and manufacturing tolerances. As it is a combination of fixed
contribution plus temperature dependent variation, most of it can be
eliminated through calibration (section 6.1).

  * •

The _cross coupling error_ or non orthogonality error is a fixed contribution
that arises from the misalignment of the sensitive axes of the inertial
sensors with respect to the orthogonal axes of the platform frame due to
manufacturing limitations, and can also be highly reduced through calibration.
The scale factor and cross coupling errors errors are in the order of
$\displaystyle\mathrm{10^{-4}}$ and
$\displaystyle\mathrm{10^{-3}\left[-\right]}$ for most inertial sensors,
although they can be higher for some low grade gyroscopes. The cross coupling
error is equal to the sine of the misalignment, which is listed by some
manufacturers.

  * •

_System noise_ or random noise is inherent to all inertial sensors and can
combine electrical, mechanical, resonance, and quantization sources. It can
originate at the sensor itself or at any other electronic equipment that
interferes with it. System noise is a stochastic process usually modeled as
white noise because its noise spectrum is approximately white, and can not be
calibrated as there is no correlation between past and future values. A white
noise process is characterized by its spectral power density PSD, which is
constant as it does not depend on the signal frequency. It corresponds to the
1 [$\displaystyle\mathrm{sec}$] crossing of the sensor Allan curve [11, 13].

System noise is sometimes referred to as random walk, which can generate
confusion with the bias. The reason is that the inertial sensor outputs are
always integrated to obtain ground velocity in the case of accelerometers and
aircraft attitude in the case of gyroscopes. As the integration of a white
noise process is indeed a random walk, the later term is commonly employed to
refer to system noise.

Table • ‣ 3.1 contains typical values for accelerometer and gyroscope root PSD
according to sensor grade [5]:

IMU Grade |  Accelerometer Root PSD $\displaystyle\mathrm{\left[m/sec/hr^{0.5}\right]}$ |  Gyroscope Root PSD $\displaystyle\mathrm{\left[deg/hr^{0.5}\right]}$
---|---|---
Aviation | 0.012 | 0.002
Tactical | 0.06 | 0.03 - 0.1
Automotive | 0.6 | 1
Table 5: Typical inertial sensor system noise according to IMU grade

  * •

Other minor error sources not considered in this article are g-dependent bias
(sensitivity of spinning mass and vibratory gyroscopes to specific force),
scale factor nonlinearity, and higher order errors (spinning mass gyroscopes
and pendulous accelerometers).

###  3.2 Single Axis Inertial Sensor Error Model

Since all sensors modeled in this article are required to provide measurements
at equispaced discrete times $\displaystyle\mathrm{t_{s}=s\cdot\Delta
t_{\scriptscriptstyle SENSED}=s\cdot\Delta t}$, this section focuses on
obtaining a discrete model for the bias and white noise errors of a single
axis inertial sensor. The results obtained here will be employed in the
following sections to generate a comprehensive IMU model.

Let’s consider a sensor in which the difference between its measurement at any
given time $\displaystyle\mathrm{\widetilde{x}\left(t\right)}$ and the real
value of the physical magnitude being measured at that same time
$\displaystyle\mathrm{x\left(t\right)}$ can be represented by a zero mean
white noise Gaussian process $\displaystyle\mathrm{\eta_{v}\left(t\right)}$
with spectral density $\displaystyle\mathrm{\sigma_{v}^{2}}$:

$\mathrm{\widetilde{x}\left(t\right)=x\left(t\right)+\eta_{v}\left(t\right)}$
(2)

Dividing (2) by $\displaystyle\mathrm{\Delta t_{\scriptscriptstyle
SENSED}=\Delta t}$ and integrating results in:

$\mathrm{\dfrac{1}{\Delta t}\int_{t_{0}}^{t_{0}+\Delta
t}\widetilde{x}\left(t\right)\,dt=\dfrac{1}{\Delta
t}\int_{t_{0}}^{t_{0}+\Delta
t}\left[x\left(t\right)+\eta_{v}\left(t\right)\right]\,dt}$ (3)

Assuming that the measurement and real value are both constant over the
integration interval121212Note that the stochastic process
$\displaystyle\mathrm{\eta_{v}}$ can not be considered constant over any
interval. [14] yields

$\mathrm{\widetilde{x}\left(t_{0}+\Delta t\right)=x\left(t_{0}+\Delta
t\right)+\dfrac{1}{\Delta t}\int_{t_{0}}^{t_{0}+\Delta
t}\eta_{v}\left(t\right)\,dt}$ (4)

This expression results in the white noise sensor error
$\displaystyle\mathrm{w\left(t\right)}$, which is the difference between the
sensor measurement $\displaystyle\mathrm{\widetilde{x}\left(t\right)}$ and the
true value $\displaystyle\mathrm{x\left(t\right)}$. Its mean and variance can
be readily computed:

$\displaystyle\displaystyle\mathrm{w\left(t_{0}+\Delta t\right)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\dfrac{1}{\Delta
t}\int_{t_{0}}^{t_{0}+\Delta t}\eta_{v}\left(t\right)\,dt}$ (5)
$\displaystyle\displaystyle\mathrm{E\left[w\left(t_{0}+\Delta
t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{0}$ (6)
$\displaystyle\displaystyle\mathrm{Var\left[w\left(t_{0}+\Delta
t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\dfrac{\sigma_{v}^{2}}{\Delta t}}$ (7)

Based on these results, the white noise error can be modeled by a discrete
random variable identically distributed to the above continuous white noise
error, this is, one that results in the same mean and variance, where
$\displaystyle\mathrm{N_{v}\sim N\left(0,\,1\right)}$ is a standard normal
random variable:

$\mathrm{w\left(s\,\Delta t_{\scriptscriptstyle
SENSED}\right)=w\left(s\,\Delta t\right)=\dfrac{\sigma_{v}}{\Delta
t^{1/2}}\,N_{vs}}$ (8)

Let’s now consider a second model in which the measurement error or bias is
given by a first order random walk process or integration of a zero mean white
noise Gaussian process $\displaystyle\mathrm{\eta_{u}\left(t\right)}$ with
spectral density $\displaystyle\mathrm{\sigma_{u}^{2}}$:

$\mathrm{\dot{b}\left(t\right)=\eta_{u}\left(t\right)\ \longrightarrow\
b\left(t_{0}+\Delta t\right)=b\left(t_{0}\right)+\int_{t_{0}}^{t_{0}+\Delta
t}\eta_{u}\left(t\right)\,dt}$ (9)

Its mean and variance can be quickly computed:

$\displaystyle\displaystyle\mathrm{E\left[b\left(t_{0}+\Delta
t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{E\left[b\left(t_{0}\right)\right]}$ (10)
$\displaystyle\displaystyle\mathrm{Var\left[b\left(t_{0}+\Delta
t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\sigma_{u}^{2}\,\Delta t}$ (11)

These results indicate that the bias can be modeled by a discrete random
variable identically distributed (this is, resulting in the same expected
value and variance) to the continuous random walk above:

$\mathrm{b\left(t_{0}+\Delta t\right)=b\left(t_{0}\right)+\sigma_{u}\,\Delta
t^{1/2}\,N_{u}}$ (12)

where $\displaystyle\mathrm{N_{u}\sim N\left(0,\,1\right)}$ is a standard
normal random variable. Operating with the above expression results in the
final expression for the discrete bias, as well as its mean and variance:

$\displaystyle\displaystyle\mathrm{b\left(s\,\Delta t\right)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{B_{0}\,N_{u0}+\sigma_{u}\,\Delta
t^{1/2}\,\sum_{i=1}^{s}N_{ui}}$ (13)
$\displaystyle\displaystyle\mathrm{E\left[b\left(s\,\Delta t\right)\right]}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{0}$ (14)
$\displaystyle\displaystyle\mathrm{Var\left[b\left(s\,\Delta t\right)\right]}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{B_{0}^{2}+\sigma_{u}^{2}\,s\,\Delta t}$
(15)

A comprehensive single axis sensor error model without scale factor can hence
be constructed by adding together the influence of the system noise provided
by (8) and the bias given by (13) [15], while assuming that the standard
normal random variables $\displaystyle\mathrm{N_{u}}$ and
$\displaystyle\mathrm{N_{v}}$ are uncorrelated131313Note that the expected
value and variance of each of the two discrete components of this sensor model
coincide with those of their continuous counterparts, but their combined mean
and variance provided by expressions (17) and (18) differ from that of the
combination of the two continuous error models given by (5) and (9). This is
the case even if considering that the two zero mean white noise Gaussian
processes $\displaystyle\mathrm{\eta_{u}}$ and
$\displaystyle\mathrm{\eta_{v}}$ are independent and hence uncorrelated. It is
however possible to obtain a discrete model whose discrete bias and white
noise components are not only identically distributed that those of their
continuous counterparts [14], even adding the equivalence of covariance
between the bias and the sensor error, but this results in a significantly
more complex model that behaves similarly to the one above at all but the
shortest time samples after sensor initialization. The author has decided not
to do so in the model described in this article, reducing complexity with
little or no loss of realism.:

$\displaystyle\displaystyle\mathrm{e_{\scriptscriptstyle BW}\left(s\,\Delta
t\right)}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\widetilde{x}\left(s\,\Delta
t\right)-x\left(s\,\Delta t\right)=B_{0}\,N_{u0}+\sigma_{u}\,\Delta
t^{1/2}\,\sum_{i=1}^{s}N_{ui}+\dfrac{\sigma_{v}}{\Delta t^{1/2}}\,N_{vs}}$
(16) $\displaystyle\displaystyle\mathrm{E\left[e_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{0}$ (17)
$\displaystyle\displaystyle\mathrm{Var\left[e_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{B_{0}^{2}+\sigma_{u}^{2}\,s\,\Delta
t+\dfrac{\sigma_{v}^{2}}{\Delta t}}$ (18)

The _discrete sensor error_ or difference between the measurement provided by
the sensor $\displaystyle\mathrm{\widetilde{x}\left(s\,\Delta
t_{\scriptscriptstyle SENSED}\right)=\widetilde{x}\left(s\,\Delta t\right)}$
at any given discrete time $\displaystyle\mathrm{s\,\Delta t}$ and the real
value of the physical variable being measured at that same discrete time
$\displaystyle\mathrm{x\left(s\,\Delta t\right)}$ is the combination of a bias
or first order random walk and a white noise process, and depends on three
parameters: the bias offset $\displaystyle\mathrm{B_{0}}$, the bias
instability $\displaystyle\mathrm{\sigma_{u}}$, and the white noise
$\displaystyle\mathrm{\sigma_{v}}$. The contributions of these three different
sources to the sensor error as well as to its first and second
integrals141414Gyroscopes measure angular velocity and their output needs to
be integrated once to obtain attitude, while accelerometers measure specific
force and are integrated once to obtain velocity and twice to obtain position.
are very different and inherent to many of the challenges encountered when
employing accelerometers and gyroscopes for inertial navigation, as explained
below.

$\displaystyle 0$$\displaystyle 200$$\displaystyle 400$$\displaystyle
600$$\displaystyle 800$$\displaystyle
1{,}000$$\displaystyle-2$$\displaystyle-1$$\displaystyle 0$$\displaystyle
1$$\displaystyle
2$$\displaystyle\mathrm{t\left[sec\right]}$$\displaystyle\mathrm{E\left[e_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]\,\left[10^{-1}\cdot unit\right]}$average 50
runstheory

$\displaystyle 0$$\displaystyle 200$$\displaystyle 400$$\displaystyle
600$$\displaystyle 800$$\displaystyle 1{,}000$$\displaystyle 0$$\displaystyle
0.2$$\displaystyle 0.4$$\displaystyle 0.6$$\displaystyle 0.8$$\displaystyle
1$$\displaystyle 1.2$$\displaystyle 1.4$$\displaystyle
1.6$$\displaystyle\mathrm{t\left[sec\right]}$$\displaystyle\mathrm{Var\left[e_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]^{1/2}\,\left[10^{-1}\cdot
unit\right]}$average 50 runstheory: totaltheory: bias offsettheory: bias
instabilitytheory: white noise

Figure 2: Propagation with time of sensor error mean and standard deviation

Figure 2 represents the performance of a fictitious sensor of
$\displaystyle\mathrm{B_{0}=1.6\cdot 10^{-2}}$,
$\displaystyle\mathrm{\sigma_{u}=4\cdot 10^{-3}}$, and
$\displaystyle\mathrm{\sigma_{v}=1\cdot 10^{-3}}$ working at a frequency of
$\displaystyle\mathrm{100\left[hz\right]\ \left(\Delta
t=0.01\left[sec\right]\right)}$, and is intended to showcase the different
behavior and relative influence on the total error of each of its three
components. The figures show the theoretical variation with time of the sensor
error mean (left side) and standard deviation (right side) given by (17) and
(18) together with the average of fifty different runs. In addition, the left
figure also includes ten of those runs to showcase the variability in results
implicit to the random variables151515Although the data is generated at
$\displaystyle\mathrm{100\left[hz\right]}$, the figure only employs 1 out of
every 1000 points, so it appears far less noisy than the real data., while the
right side shows the theoretical contribution to the standard deviation of
each of the three components. In addition to the near equivalence between the
theory and the average of several runs, the figure shows that the bias
instability is the commanding long term factor in the deviation between the
sensor measurement and its zero mean (the standard deviation of the bias
instability grows with the square root of time while the other two components
are constant). As discussed in section 3.1, the bias drift or bias instability
is indeed the most important quality parameter of an inertial sensor. This is
also the case when the sensor output is integrated, as discussed below.

Let’s integrate the sensor error over a timespan
$\displaystyle\mathrm{s\,\Delta t}$ to evaluate the growth with time of both
its expected value and its variance161616As the interest lies primarily in
$\displaystyle\mathrm{s\gg 1}$, a simple integration method such as the
rectangular rule is employed.:

$\displaystyle\displaystyle\mathrm{f_{\scriptscriptstyle BW}\left(s\,\Delta
t\right)}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{f_{\scriptscriptstyle
BW}\left(0\right)+\int_{0}^{s\,\Delta t}e_{\scriptscriptstyle
BW}\left(\tau\right)\,d\tau=f_{\scriptscriptstyle BW}\left(0\right)+\Delta
t\sum_{i=1}^{s}e_{\scriptscriptstyle BW}\left(i\,\Delta t\right)}$ (19)
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{f_{\scriptscriptstyle
BW}\left(0\right)+B_{0}\,N_{u0}\,s\,\Delta t+\sigma_{u}\,\Delta
t^{3/2}\sum_{i=1}^{s}\left(s-i+1\right)\,N_{ui}+\sigma_{v}\,\Delta
t^{1/2}\sum_{i=1}^{s}N_{vi}}$
$\displaystyle\displaystyle\mathrm{E\left[f_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{f_{\scriptscriptstyle BW}\left(0\right)}$
(20) $\displaystyle\displaystyle\mathrm{Var\left[f_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{B_{0}^{2}\,\left(s\,\Delta
t\right)^{2}+\dfrac{\sigma_{u}^{2}}{6}\,\Delta
t^{3}\,s\,\left(s+1\right)\,\left(2\,s+1\right)+\sigma_{v}^{2}\,\left(s\,\Delta
t\right)}$ (21) $\displaystyle\displaystyle\mathrm{\approx}$
$\displaystyle\displaystyle\mathrm{B_{0}^{2}\,\left(s\,\Delta
t\right)^{2}+\dfrac{\sigma_{u}^{2}}{3}\,\left(s\,\Delta
t\right)^{3}+\sigma_{v}^{2}\,\left(s\,\Delta t\right)}$

$\displaystyle 0$$\displaystyle 200$$\displaystyle 400$$\displaystyle
600$$\displaystyle 800$$\displaystyle
1{,}000$$\displaystyle-100$$\displaystyle-50$$\displaystyle 0$$\displaystyle
50$$\displaystyle 100$$\displaystyle
150$$\displaystyle\mathrm{t\left[sec\right]}$$\displaystyle\mathrm{E\left[f_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]\,\left[unit\cdot sec\right]}$average 50
runstheory

$\displaystyle 0$$\displaystyle 200$$\displaystyle 400$$\displaystyle
600$$\displaystyle 800$$\displaystyle 1{,}000$$\displaystyle 0$$\displaystyle
20$$\displaystyle 40$$\displaystyle 60$$\displaystyle
80$$\displaystyle\mathrm{t\left[sec\right]}$$\displaystyle\mathrm{Var\left[f_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]^{1/2}\,\left[unit\cdot sec\right]}$average
50 runstheory: totaltheory: bias offsettheory: bias instabilitytheory: white
noise

Figure 3: Propagation with time of first integral of sensor error mean and
standard deviation

Figure 3 follows the same pattern as figure 2 but applied to the error
integral instead of to the error itself. It would represent the attitude error
resulting from integrating the gyroscope output, or the velocity error
expected when integrating the specific force measured by an accelerometer. The
conclusions are the same as before but significantly more accentuated. Not
only is the expected value of the error constant instead of zero
($\displaystyle\mathrm{f_{\scriptscriptstyle BW}\left(0\right)=3}$ has been
employed in the experiment), but the growth in the standard deviation (over a
non zero mean) is much quicker than before. The bias instability continues to
be the dominating factor but now increases with a power of
$\displaystyle\mathrm{t^{3/2}}$, while the bias offset and white noise
contributions also increase with time, although with powers of
$\displaystyle\mathrm{t}$ and $\displaystyle\mathrm{t^{1/2}}$ respectively.
Let’s continue the process and integrate the error a second time:

$\displaystyle\displaystyle\mathrm{g_{\scriptscriptstyle BW}\left(s\,\Delta
t\right)}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{g_{\scriptscriptstyle
BW}\left(0\right)+\int_{0}^{s\,\Delta t}f_{\scriptscriptstyle
BW}\left(\tau\right)\,d\tau=g_{\scriptscriptstyle BW}\left(0\right)+\Delta
t\sum_{i=1}^{s}f_{\scriptscriptstyle BW}\left(i\,\Delta
t\right)=g_{\scriptscriptstyle BW}\left(0\right)+f_{\scriptscriptstyle
BW}\left(0\right)\,s\,\Delta t+}$ (22)
$\displaystyle\displaystyle\mathrm{+\dfrac{B_{0}}{2}\,N_{u0}\,\left(s\,\Delta
t\right)^{2}+\sigma_{u}\,\Delta
t^{5/2}\sum_{i=1}^{s}\sum_{j=1}^{s-i+1}j\,N_{ui}+\sigma_{v}\,\Delta
t^{3/2}\sum_{i=1}{s}\left(s-i+1\right)\,N_{vi}}$
$\displaystyle\displaystyle\mathrm{E\left[g_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{g_{\scriptscriptstyle
BW}\left(0\right)+f_{\scriptscriptstyle BW}\left(0\right)\,\left(s\,\Delta
t\right)}$ (23)
$\displaystyle\displaystyle\mathrm{Var\left[g_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\dfrac{B_{0}^{2}}{4}\,\left(s\,\Delta
t\right)^{4}+\sigma_{u}^{2}\,\Delta
t^{5}\,\sum_{i=1}^{s}\left(\sum_{j=1}^{s-i+1}j\right)^{2}+\dfrac{\sigma_{v}^{2}}{6}\,\Delta
t^{3}\,s\,\left(s+1\right)\,\left(2\,s+1\right)}$ (24)
$\displaystyle\displaystyle\mathrm{\approx}$
$\displaystyle\displaystyle\mathrm{\dfrac{B_{0}^{2}}{4}\,\left(s\,\Delta
t\right)^{4}+\dfrac{\sigma_{u}^{2}}{20}\,\left(s\,\Delta
t\right)^{5}+\dfrac{\sigma_{v}^{2}}{3}\,\left(s\,\Delta t\right)^{3}}$

Figure 4 shows the same type of figures but applied to the second integral of
the error ($\displaystyle\mathrm{g_{\scriptscriptstyle BW}\left(0\right)=1.5}$
has been employed in the experiment). In this case the degradation of the
results with time is even more intense to the point were the measurements are
useless after a very short period of time. Unless corrected by the navigation
system, this is equivalent to the error in position obtained by double
integrating the output of the accelerometers.

$\displaystyle 0$$\displaystyle 200$$\displaystyle 400$$\displaystyle
600$$\displaystyle 800$$\displaystyle
1{,}000$$\displaystyle-4$$\displaystyle-2$$\displaystyle 0$$\displaystyle
2$$\displaystyle 4$$\displaystyle
6$$\displaystyle\mathrm{t\left[sec\right]}$$\displaystyle\mathrm{E\left[g_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]\,\left[10^{4}\cdot unit\cdot
sec^{2}\right]}$average 50 runstheory

$\displaystyle 0$$\displaystyle 200$$\displaystyle 400$$\displaystyle
600$$\displaystyle 800$$\displaystyle 1{,}000$$\displaystyle 0$$\displaystyle
0.5$$\displaystyle 1$$\displaystyle 1.5$$\displaystyle 2$$\displaystyle
2.5$$\displaystyle
3$$\displaystyle\mathrm{t\left[sec\right]}$$\displaystyle\mathrm{Var\left[g_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)\right]^{1/2}\,\left[10^{4}\cdot unit\cdot
sec^{2}\right]}$average 50 runstheory: totaltheory: bias offsettheory: bias
instabilitytheory: white noise

Figure 4: Propagation with time of second integral of sensor error mean and
standard deviation

Let’s summarize the main points of the single axis inertial sensor discrete
error model developed in this section, that includes the influence of the bias
and the system error but not that of the scale factor and cross coupling
errors included in the three dimensional error model of section 3.9. The error
$\displaystyle\mathrm{e_{\scriptscriptstyle BW}\left(s\,\Delta t\right)}$,
which applies to specific force for accelerometers and inertial angular
velocity in the case of gyroscopes, depends on three factors: bias offset
$\displaystyle\mathrm{B_{0}}$, bias drift $\displaystyle\mathrm{\sigma_{u}}$,
and white noise $\displaystyle\mathrm{\sigma_{v}}$. Its mean is always zero,
but the error standard deviation grows with time
($\displaystyle\mathrm{\propto\,t^{1/2}}$) due to the bias drift with constant
contributions from the bias offset and the white noise. When integrating the
error to obtain $\displaystyle\mathrm{f_{\scriptscriptstyle BW}\left(s\,\Delta
t\right)}$, equivalent to ground velocity for accelerometers and attitude for
gyroscopes, the initial speed error or initial attitude error
$\displaystyle\mathrm{f_{\scriptscriptstyle BW}\left(0\right)}$ becomes the
fourth contributor, and an important one indeed, as it becomes the mean of the
first integral error at any time. The standard deviation, which measures the
spread over the non zero mean, increases very quickly with time because of the
bias instability ($\displaystyle\mathrm{\propto\,t^{3/2}}$), with
contributions from the offset ($\displaystyle\mathrm{\propto\,t}$) and the
white noise ($\displaystyle\mathrm{\propto\,t^{1/2}}$). If integrating a
second time to obtain $\displaystyle\mathrm{g_{\scriptscriptstyle
BW}\left(s\,\Delta t\right)}$, or position in case of the accelerometer, the
initial position error $\displaystyle\mathrm{g_{\scriptscriptstyle
BW}\left(0\right)}$ turns into the fifth contributor. The expected value of
the position error grows linearly with time due to the initial velocity error
with an additional constant contribution from the initial position error,
while the position standard deviation (measuring spread over a growing average
value) grows extremely quick due mostly to the bias instability
($\displaystyle\mathrm{\propto\,t^{5/2}}$), but also because of the bias
offset ($\displaystyle\mathrm{\propto\,t^{2}}$) and the white noise
($\displaystyle\mathrm{\propto\,t^{3/2}}$). Table 3.2 shows the standard units
of the different sources of error for both accelerometers and gyroscopes.

Units | $\displaystyle\mathrm{B_{0}}$ | $\displaystyle\mathrm{\sigma_{u}}$ | $\displaystyle\mathrm{\sigma_{v}}$ | $\displaystyle\mathrm{f_{\scriptscriptstyle BW}\left(0\right)}$ | $\displaystyle\mathrm{g_{\scriptscriptstyle BW}\left(0\right)}$
---|---|---|---|---|---
Accelerometer | $\displaystyle\mathrm{m/sec^{2}}$ | $\displaystyle\mathrm{m/sec^{2.5}}$ | $\displaystyle\mathrm{m/sec^{1.5}}$ | $\displaystyle\mathrm{m/sec}$ | $\displaystyle\mathrm{m}$
Gyroscope | $\displaystyle\mathrm{deg/sec}$ | $\displaystyle\mathrm{deg/sec^{1.5}}$ | $\displaystyle\mathrm{deg/sec^{0.5}}$ | $\displaystyle\mathrm{deg}$ | N/A
Table 6: Units for single axis inertial sensor error sources

###  3.3 Obtainment of System Noise Values

This section focuses on the significance of system or white noise error
$\displaystyle\mathrm{\sigma_{v}}$ and how to obtain it from sensor
specifications, which often refer to the integral of the output instead of the
output itself. As the integral of white noise is a random walk process, the
angle random walk of a gyroscope is equivalent to white noise in the angular
rate output, while velocity random walk refers to specific force white noise
in accelerometers [8]. The discussion that follows applies to gyroscopes but
is fully applicable to accelerometers if replacing angular rate by specific
force and attitude or angle by ground velocity.

Angle random walk, measured in [$\displaystyle\mathrm{rad/sec^{1/2}}$],
[$\displaystyle\mathrm{deg/hr^{1/2}}$], or equivalent, describes the average
deviation or error that occurs when the sensor output signal is integrated due
to system noise exclusively, without considering other error sources such as
bias or scale factor [16]. If integrating multiple times and obtaining a
distribution of end points at a given final time
$\displaystyle\mathrm{s\,\Delta t}$, the standard deviation of this
distribution, containing the final angles at the final time, scales linearly
with the white noise level $\displaystyle\mathrm{\sigma_{v}}$, the square root
of the integration step size $\displaystyle\mathrm{\Delta t}$, and the square
root of the number of steps $\displaystyle\mathrm{s}$, as noted by the last
term of (21). This means that an angle random walk of 1
[$\displaystyle\mathrm{deg/sec^{1/2}}$] translates into an standard deviation
for the error of 1 [$\displaystyle\mathrm{deg}$] after 1
[$\displaystyle\mathrm{sec}$], 10 [$\displaystyle\mathrm{deg}$] after 100
[$\displaystyle\mathrm{sec}$], and $\displaystyle\mathrm{1000^{1/2}\approx
31.6}$ [$\displaystyle\mathrm{deg}$] after 1000 [$\displaystyle\mathrm{sec}$].

Manufacturers often provide this information as the power spectral density PSD
of the white noise process in [$\displaystyle\mathrm{deg^{2}/hr^{2}/hz}$] or
equivalent, where it is necessary to take its square root to obtain
$\displaystyle\mathrm{\sigma_{v}}$, or as the root PSD in
[$\displaystyle\mathrm{deg/hr/hz^{1/2}}$] that is equivalent to
$\displaystyle\mathrm{\sigma_{v}}$. Sometimes it is even provided as the PSD
of the random walk process, not the white noise, in units
[$\displaystyle\mathrm{deg/hr}$] or equivalent. It is then necessary to
multiply this number by the square root of the sampling interval
$\displaystyle\mathrm{\Delta t}$ or divide it by the square root of the
sampling frequency to obtain the desired $\displaystyle\mathrm{\sigma_{v}}$
value.

###  3.4 Obtainment of Bias Drift Values

This section describes the meaning of bias instability
$\displaystyle\mathrm{\sigma_{u}}$ (also known as bias stability or bias
drift) and how to obtain it from sensor specifications. As in the previous
section, the discussion is centered on gyroscopes but it is fully applicable
to accelerometers as well. Bias instability can be defined as the potential of
the sensor error to stay within a certain range for a certain time [17]. A
small number of manufacturers directly provide sensor output changes over
time, which directly relates with the bias instability (also known as in-run
bias variation, bias drift, or rate random walk) per the second term of (18).
If provided with an angular rate change of x [$\displaystyle\mathrm{deg/sec}$]
($\displaystyle\mathrm{1\sigma}$) in t [$\displaystyle\mathrm{sec}$], then
$\displaystyle\mathrm{\sigma_{u}}$ can be obtained as follows [14, 18]:

$\mathrm{\sigma_{u}=\dfrac{x}{t^{1/2}}}$ (25)

As the bias drift is responsible for the growth of sensor error with time
(figure 2), manufacturers more commonly provide bias stability measurements
that describe how the bias of a device may change over a specified period of
time [15], typically around 100 [$\displaystyle\mathrm{sec}$]. Bias stability
is usually specified as a $\displaystyle\mathrm{1\sigma}$ value with units
[$\displaystyle\mathrm{deg/hr}$] or [$\displaystyle\mathrm{deg/sec}$], which
can be interpreted as follows according to (16, 17, 18). If the sensor error
(or bias) is known at a given time t, then a $\displaystyle\mathrm{1\sigma}$
bias stability of 0.01 [$\displaystyle\mathrm{deg/hr}$] over 100
[$\displaystyle\mathrm{sec}$] means that the bias at time
$\displaystyle\mathrm{t+100}$ [$\displaystyle\mathrm{sec}$] is a random
variable with mean the error at time t and standard deviation 0.01
[$\displaystyle\mathrm{deg/hr}$], and expression (25) can be used to obtain
$\displaystyle\mathrm{\sigma_{u}}$. As the bias behaves as a random walk over
time whose standard deviation grows proportionally to the square root of time,
the bias stability is sometimes referred as a bias random walk.

In reality bias fluctuations do not really behave as a random walk. If they
did, the uncertainty in the bias of a device would grow without bound as the
timespan increased, which is not the case. In practice the bias is constrained
to be within some range, and therefore the random walk model is only a good
approximation to the true process for short periods of time [15].

###  3.5 Platform, Accelerometers, and Gyroscopes Frames

The following sections make use of three different reference frames to
describe the readings of accelerometers and gyroscopes:

  * •

The _platform frame_ $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ is a
Cartesian reference system with its origin located at the IMU reference point
(section 3.8) and its three axes
$\displaystyle\mathrm{\left\\{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle
P},\,\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle
P},\,\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle P}\right\\}}$ forming a right hand system that is
loosely aligned with the aircraft body axes, so they point in the general
directions of the aircraft fuselage (forward), aircraft wings (rightwards),
and downward, respectively [3, 9, 19].

A proper definition of the platform frame is indispensable for navigation, as
the calibrated outputs of the accelerometers and gyroscopes are based on it
(sections 3.6 and 3.7). The $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$
frame can be obtained from the body frame
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ by a rotation best described
by the Euler angles
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}^{\scriptscriptstyle
BP}=\left[\psi^{\scriptscriptstyle P},\ \theta^{\scriptscriptstyle P},\
\xi^{\scriptscriptstyle P}\right]^{T}}$171717These Euler angles correspond to
the 3 - 2 - 1 (yaw, pitch, roll) convention employed in aeronautics. followed
by a translation
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}$ (section 3.8) from the aircraft center of mass to
the IMU reference point.

  * •

The _accelerometers frame_ $\displaystyle\mathrm{F_{\scriptscriptstyle A}}$ is
a non orthogonal reference system also centered at the IMU reference point [3,
9, 19]. The basis vectors
$\displaystyle\mathrm{\left\\{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle
A},\,\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle
A},\,\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle A}\right\\}}$ are aligned with each of the three
accelerometers sensing axes181818Each accelerometer hence only senses the
specific force component parallel to its sensing axis. (section 3.6), but they
are not orthogonal among them due to manufacturing inaccuracies. This implies
that the angles between the $\displaystyle\mathrm{F_{\scriptscriptstyle A}}$
and $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ axes are very small.

It is always possible, with no loss of generality, to impose that
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle P}}$ coincides with
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle A}}$ and that
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle P}}$ is located in the plane defined by
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle A}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle A}}$. If this is the case,
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle
A}\perp\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle P}}$,
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle
A}\perp\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle P}}$, and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle
A}\perp\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle P}}$, and the relative attitude between the
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ and
$\displaystyle\mathrm{F_{\scriptscriptstyle A}}$ frames can be defined by
three independent small rotations.

    * –

The
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle A}}$ axis can be obtained from
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle P}}$ by means of a small rotation
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle{ACC,3}}}$ about
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle P}}$.

    * –

The
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle A}}$ axis can be obtained from
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle P}}$ by two small rotations:
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle{ACC,1}}}$ about
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle P}}$ and
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle{ACC,2}}}$ about
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle P}}$.

Although the exact relationships can be obtained [9], and given that the
angles are very small, it is possible to consider
$\displaystyle\mathrm{\cos\alpha_{\scriptscriptstyle{ACC,i}}=1,\
\sin\alpha_{\scriptscriptstyle{ACC,i}}=\alpha_{\scriptscriptstyle{ACC,i}}}$,
and
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle{ACC,i}}\cdot\alpha_{\scriptscriptstyle{ACC,j}}=0\
\forall\ i,j\in\left\\{1,\ 2,\ 3\right\\},i\neq j}$, resulting in the
following transformations between free vectors viewed in the platform
($\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
P}}$) and accelerometer
($\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
A}}$) frames respectively191919As $\displaystyle\mathrm{F_{\scriptscriptstyle
A}}$ is not orthogonal, the transformation matrices are denoted with
$\displaystyle\mathrm{\star}$ to indicate that they are not proper rotation
matrices.:

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
P}=\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PA}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
A}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\begin{bmatrix}1&0&0\\\
\mathrm{\alpha_{\scriptscriptstyle{ACC,3}}}&1&0\\\
-\mathrm{\alpha_{\scriptscriptstyle{ACC,2}}}&\mathrm{\alpha_{\scriptscriptstyle{ACC,1}}}&1\end{bmatrix}\
\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle A}}$ (26)
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
A}=\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
AP}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
P}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\begin{bmatrix}1&0&0\\\
-\mathrm{\alpha_{\scriptscriptstyle{ACC,3}}}&1&0\\\
\mathrm{\alpha_{\scriptscriptstyle{ACC,2}}}&-\mathrm{\alpha_{\scriptscriptstyle{ACC,1}}}&1\end{bmatrix}\
\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle P}}$ (27)

  * •

The _gyroscopes frame_ $\displaystyle\mathrm{F_{\scriptscriptstyle Y}}$ is
similar to the accelerometers frame
$\displaystyle\mathrm{F_{\scriptscriptstyle A}}$ defined above, but aligned
with the gyroscopes sensing axes instead of those of the accelerometers [3, 9,
19]. It is also a non orthogonal reference system centered at the IMU
reference point, but no simplifications can be made about the relative
orientation of its axes
$\displaystyle\mathrm{\left\\{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
1}^{\scriptscriptstyle
Y},\,\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
2}^{\scriptscriptstyle
Y},\,\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{\scriptscriptstyle
3}^{\scriptscriptstyle Y}\right\\}}$ with respect to those of
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$, so their relative attitude
is defined by six small rotations
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle{GYR,ij}}\ \forall\
i,j\in\left\\{1,\ 2,\ 3\right\\},i\neq j}$, where
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle{GYR,ij}}}$ is the rotation of
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{i}^{\scriptscriptstyle
Y}}$ about
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{i}$}}_{j}^{\scriptscriptstyle
P}}$.

An approach similar to that employed with accelerometers leads to the
following transformations between free vectors viewed in the platform
($\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
P}}$) and gyroscope
($\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
Y}}$) frames:

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
P}=\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PY}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
Y}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\begin{bmatrix}1&-\mathrm{\alpha_{\scriptscriptstyle{GYR,23}}}&\mathrm{\alpha_{\scriptscriptstyle{GYR,32}}}\\\
\mathrm{\alpha_{\scriptscriptstyle{GYR,13}}}&1&-\mathrm{\alpha_{\scriptscriptstyle{GYR,31}}}\\\
-\mathrm{\alpha_{\scriptscriptstyle{GYR,12}}}&\mathrm{\alpha_{\scriptscriptstyle{GYR,21}}}&1\end{bmatrix}\
\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle Y}}$ (28)
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
Y}=\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
YP}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle
P}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\begin{bmatrix}1&\mathrm{\alpha_{\scriptscriptstyle{GYR,23}}}&-\mathrm{\alpha_{\scriptscriptstyle{GYR,32}}}\\\
-\mathrm{\alpha_{\scriptscriptstyle{GYR,13}}}&1&\mathrm{\alpha_{\scriptscriptstyle{GYR,31}}}\\\
\mathrm{\alpha_{\scriptscriptstyle{GYR,12}}}&-\mathrm{\alpha_{\scriptscriptstyle{GYR,21}}}&1\end{bmatrix}\
\mbox{\boldmath{$\displaystyle\mathrm{v}$}}^{\scriptscriptstyle P}}$ (29)

###  3.6 Accelerometer Triad Sensor Error Model

An IMU is equipped with an accelerometer triad composed by three individual
accelerometers, each of which measures the projection of the specific force
over its sensing axis as described in section 3.2 while incurring in an error
$\displaystyle\mathrm{e_{\scriptscriptstyle BW,ACC}}$ that can be modeled as a
combination of bias offset, bias drift, and white noise (16). The three
accelerometers can be considered infinitesimally small and located at the _IMU
reference point_ , which is defined as the intersection between the sensing
axes of the three accelerometers. As the accelerometer frame
$\displaystyle\mathrm{F_{\scriptscriptstyle A}}$ is centered at the IMU
reference point and its three non orthogonal axes coincide with the
accelerometers sensing axes, (30) joins together the measurements of the three
individual accelerometers:

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IA}^{\scriptscriptstyle
A}=\mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle ACC}\
\left(\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IA}^{\scriptscriptstyle
A}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle A}\right)}$ (30)

where
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IA}^{\scriptscriptstyle A}}$ is the specific force viewed in the accelerometer
frame $\displaystyle\mathrm{F_{\scriptscriptstyle A}}$,
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IA}^{\scriptscriptstyle A}}$ represents its measurement also viewed in
$\displaystyle\mathrm{F_{\scriptscriptstyle A}}$,
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle A}}$ is the error introduced by each accelerometer
(16), and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle
ACC}}$ is a square diagonal matrix containing the scale factor errors
$\displaystyle\mathrm{\left\\{s_{\scriptscriptstyle{ACC,1}},\
s_{\scriptscriptstyle{ACC,2}},\ s_{\scriptscriptstyle{ACC,3}}\right\\}}$ for
each accelerometer (section 3.1). It is however preferred to obtain an
expression in which the specific forces are viewed in the orthogonal platform
frame $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ instead of the
accelerometer frame $\displaystyle\mathrm{F_{\scriptscriptstyle A}}$. As both
share the same origin,

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle
P}=\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PA}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle
ACC}\ \left(\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
AP}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle
P}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle A}\right)}$ (31)

where
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PA}^{\star}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
AP}^{\star}}$, defined by (26) and (27), contain the cross coupling errors
$\displaystyle\mathrm{\left\\{\alpha_{\scriptscriptstyle{ACC,1}},\
\alpha_{\scriptscriptstyle{ACC,2}},\
\alpha_{\scriptscriptstyle{ACC,3}}\right\\}}$ generated by the misalignment of
the accelerometer sensing axes. The scale factor and cross coupling errors
contain fixed and temperature dependent error contributions (refer to section
3.1) that can be modeled as normal random variables:

$\displaystyle\displaystyle\mathrm{s_{\scriptscriptstyle{ACC,i}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(1,\
s_{\scriptscriptstyle ACC}^{2}\right)\;\hskip 20.0pt\forall\ i\in\left\\{1,\
2,\ 3\right\\}}$ (32)
$\displaystyle\displaystyle\mathrm{\alpha_{\scriptscriptstyle{ACC,i}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(0,\
\alpha_{\scriptscriptstyle ACC}^{2}\right)\hskip 20.0pt\forall\
i\in\left\\{1,\ 2,\ 3\right\\}}$ (33)

where $\displaystyle\mathrm{s_{\scriptscriptstyle ACC}}$ and
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle ACC}}$ can be obtained from
the sensor specifications. Equation (31) can be transformed to make it more
useful by defining the accelerometer scale and cross coupling error matrix
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
ACC}}$:

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
ACC}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PA}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle
ACC}\ \mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
AP}^{\star}=\begin{bmatrix}\mathrm{m_{\scriptscriptstyle{ACC,11}}}&0&0\\\
\mathrm{m_{\scriptscriptstyle{ACC,21}}}&\mathrm{m_{\scriptscriptstyle{ACC,22}}}&0\\\
\mathrm{m_{\scriptscriptstyle{ACC,31}}}&\mathrm{m_{\scriptscriptstyle{ACC,32}}}&\mathrm{m_{\scriptscriptstyle{ACC,33}}}\end{bmatrix}}$
(34) $\displaystyle\displaystyle\mathrm{\approx}$
$\displaystyle\displaystyle\mathrm{\begin{bmatrix}\mathrm{s_{\scriptscriptstyle{ACC,1}}}&0&0\\\
\mathrm{\alpha_{\scriptscriptstyle{ACC,3}}\
\left(s_{\scriptscriptstyle{ACC,1}}-s_{\scriptscriptstyle{ACC,2}}\right)}&\mathrm{s_{\scriptscriptstyle{ACC,2}}}&0\\\
\mathrm{\alpha_{\scriptscriptstyle{ACC,2}}\
\left(s_{\scriptscriptstyle{ACC,3}}-s_{\scriptscriptstyle{ACC,1}}\right)}&\mathrm{\alpha_{\scriptscriptstyle{ACC,1}}\
\left(s_{\scriptscriptstyle{ACC,2}}-s_{\scriptscriptstyle{ACC,3}}\right)}&\mathrm{s_{\scriptscriptstyle{ACC,3}}}\end{bmatrix}}$

Considering that the scale and cross coupling errors are uncorrelated and very
small, and taking into account the expressions for the mean and variance of
the sum and product of two random variables [20], the different components
$\displaystyle\mathrm{m_{\scriptscriptstyle{ACC,ij}}}$ of
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
ACC}}$ can be obtained as follows $\displaystyle\mathrm{\forall\ i,\
j\in\left\\{1,\ 2,\ 3\right\\}}$:

$\displaystyle\displaystyle\mathrm{m_{\scriptscriptstyle{ACC,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(1,\
s_{\scriptscriptstyle ACC}^{2}\right)\hskip 146.0pti=j}$ (35)
$\displaystyle\displaystyle\mathrm{m_{\scriptscriptstyle{ACC,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(0,\
\left[\sqrt{2}\ \alpha_{\scriptscriptstyle ACC}\ s_{\scriptscriptstyle
ACC}\right]^{2}\right)=N\left(0,m_{\scriptscriptstyle ACC}^{2}\right)\hskip
20.0pti>j}$ (36)
$\displaystyle\displaystyle\mathrm{m_{\scriptscriptstyle{ACC,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{0\hskip
191.0pti<j}$ (37)

Let’s also define the accelerometer error transformation matrix
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{\scriptscriptstyle
ACC}}$ as

$\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{\scriptscriptstyle
ACC}=\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PA}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle
ACC}=\begin{bmatrix}\mathrm{n_{\scriptscriptstyle{ACC,11}}}&0&0\\\
\mathrm{n_{\scriptscriptstyle{ACC,21}}}&\mathrm{n_{\scriptscriptstyle{ACC,22}}}&0\\\
\mathrm{n_{\scriptscriptstyle{ACC,31}}}&\mathrm{n_{\scriptscriptstyle{ACC,32}}}&\mathrm{n_{\scriptscriptstyle{ACC,33}}}\end{bmatrix}=\begin{bmatrix}\mathrm{s_{\scriptscriptstyle{ACC,1}}}&0&0\\\
\mathrm{\alpha_{\scriptscriptstyle{ACC,3}}\
s_{\scriptscriptstyle{ACC,1}}}&\mathrm{s_{\scriptscriptstyle{ACC,2}}}&0\\\
-\mathrm{\alpha_{\scriptscriptstyle{ACC,2}}\
s_{\scriptscriptstyle{ACC,1}}}&\mathrm{\alpha_{\scriptscriptstyle{ACC,1}}\
s_{\scriptscriptstyle{ACC,2}}}&\mathrm{s_{\scriptscriptstyle{ACC,3}}}\end{bmatrix}}$
(38)

A process similar to that employed above leads to:

$\displaystyle\displaystyle\mathrm{n_{\scriptscriptstyle{ACC,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(1,\
s_{\scriptscriptstyle ACC}^{2}\right)\hskip 132.0pti=j}$ (39)
$\displaystyle\displaystyle\mathrm{n_{\scriptscriptstyle{ACC,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(0,\
\alpha_{\scriptscriptstyle ACC}^{2}\left(1+s_{\scriptscriptstyle
ACC}^{2}\right)\right)\approx N\left(0,\alpha_{\scriptscriptstyle
ACC}^{2}\right)\hskip 20.0pti>j}$ (40)
$\displaystyle\displaystyle\mathrm{n_{\scriptscriptstyle{ACC,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{0\hskip
178.0pti<j}$ (41)

Taking into account the expressions for the mean and variance of the sum and
product of two random variables [20], and knowing that the cross coupling
errors are very small $\displaystyle\mathrm{\left(1+\alpha_{\scriptscriptstyle
ACC}^{2}\approx 1\right)}$, it can be proven that the bias and white noise
errors viewed in the platform frame
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ respond to a expression
similar to (16):

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle P}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle P}\left(s\,\Delta t_{\scriptscriptstyle
SENSED}\right)=\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle P}\left(s\,\Delta t\right)}$ (42)
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{\scriptscriptstyle
ACC}\ \mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle
A}=B_{0\scriptscriptstyle{ACC}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{u0,\scriptscriptstyle
ACC}+\sigma_{u\scriptscriptstyle{ACC}}\,\Delta
t^{1/2}\,\sum_{i=1}^{s}\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{ui,\scriptscriptstyle
ACC}+\dfrac{\sigma_{v\scriptscriptstyle{ACC}}}{\Delta
t^{1/2}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{vs,\scriptscriptstyle
ACC}}$

where each
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{u,\scriptscriptstyle
ACC}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{v,\scriptscriptstyle
ACC}}$ is a random vector composed by three independent standard normal random
variables $\displaystyle\mathrm{N\left(0,\,1\right)}$. Note that as the bias
drift is mostly a warm up process that stabilizes itself after a few minutes
of operation, the random walk within (42) is not allowed to vary freely but is
restricted to within a band of width $\displaystyle\mathrm{\pm\
100\,\sigma_{u\scriptscriptstyle{ACC}}\,\Delta t^{1/2}}$. The final model for
the accelerometer measurements viewed in
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ results in

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle
P}=\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle ACC}\
\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle
P}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle P}}$ (43)

where
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
ACC}}$ is described in (34) through (37) and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle P}}$ is provided by (42). This model relies on
inputs for the bias offset
$\displaystyle\mathrm{B_{0\scriptscriptstyle{ACC}}}$, bias drift
$\displaystyle\mathrm{\sigma_{u\scriptscriptstyle{ACC}}}$, white noise
$\displaystyle\mathrm{\sigma_{v\scriptscriptstyle{ACC}}}$, scale factor error
$\displaystyle\mathrm{s_{\scriptscriptstyle ACC}}$, and cross coupling error
$\displaystyle\mathrm{m_{\scriptscriptstyle ACC}}$, which are obtained from
the accelerometer specifications shown in table 3.10 within section 3.10.

###  3.7 Gyroscopes Triad Sensor Error Model

The IMU is also equipped with a triad of gyroscopes, each of which measures
the projection of the inertial angular velocity over its sensing axis as
described in section 3.2. The obtainment of the gyroscope triad model is fully
analogous to that of the accelerometers in the previous section, with the
added difficulty that the transformation between the gyroscope frame
$\displaystyle\mathrm{F_{\scriptscriptstyle Y}}$ and platform frame
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ relies on six small angles
instead of three. The starting point hence is:

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle
P}=\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PY}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle
GYR}\ \left(\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
YP}^{\star}\
\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle
P}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle Y}\right)}$ (44)

where
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle P}}$ is the inertial angular velocity viewed in the
platform frame $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$,
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle P}}$ represents its measurement also viewed in
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$,
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle Y}}$ is the error introduced by each gyroscope
(16),
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle
GYR}}$ is a square diagonal matrix containing the scale factor errors
$\displaystyle\mathrm{\left\\{s_{\scriptscriptstyle{GYR,1}},\
s_{\scriptscriptstyle{GYR,2}},\ s_{\scriptscriptstyle{GYR,2}}\right\\}}$, and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PY}^{\star}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
YP}^{\star}}$, defined by (28) and (29), contain the cross coupling errors
$\displaystyle\mathrm{\alpha_{\scriptscriptstyle{GYR,12}},\
\alpha_{\scriptscriptstyle{GYR,21}},\ \alpha_{\scriptscriptstyle{GYR,13}},\
\alpha_{\scriptscriptstyle{GYR,31}},\ \alpha_{\scriptscriptstyle{GYR,23}},\
\alpha_{\scriptscriptstyle{GYR,32}}}$ generated by the misalignment of the
gyroscope sensing axes.

Operating in the same way as in section 3.6 leads to:

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle P}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle P}\left(s\,\Delta t_{\scriptscriptstyle
SENSED}\right)=\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle P}\left(s\,\Delta t\right)}$ (45)
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{B_{0\scriptscriptstyle{GYR}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{u0,\scriptscriptstyle
GYR}+\sigma_{u\scriptscriptstyle{GYR}}\,\Delta
t^{1/2}\,\sum_{i=1}^{s}\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{ui,\scriptscriptstyle
GYR}+\dfrac{\sigma_{v\scriptscriptstyle{GYR}}}{\Delta
t^{1/2}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{vs,\scriptscriptstyle
GYR}}$
$\displaystyle\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle P}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
GYR}\ \mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle
P}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle P}}$ (46)

where each
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{ui,\scriptscriptstyle
GYR}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{v,\scriptscriptstyle
GYR}}$ is a random vector composed by three independent standard normal random
variables $\displaystyle\mathrm{N\left(0,\,1\right)}$. As in the case of the
accelerometers, the random walk within (45) representing the bias drift is not
allowed to vary freely but is restricted to within a band of width
$\displaystyle\mathrm{\pm\ 100\,\sigma_{u\scriptscriptstyle{GYR}}\,\Delta
t^{1/2}}$. This model relies on inputs for the bias offset
$\displaystyle\mathrm{B_{0\scriptscriptstyle{GYR}}}$, bias drift
$\displaystyle\mathrm{\sigma_{u\scriptscriptstyle{GYR}}}$, white noise
$\displaystyle\mathrm{\sigma_{v\scriptscriptstyle{GYR}}}$, scale factor error
$\displaystyle\mathrm{s_{\scriptscriptstyle GYR}}$, and cross coupling error
$\displaystyle\mathrm{m_{\scriptscriptstyle GYR}}$, which are obtained from
the gyroscope specifications shown in table 3.10 within section 3.10. The
gyroscope scale and cross coupling error matrix
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
GYR}}$ responds to:

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
GYR}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
PY}^{\star}\ \mbox{\boldmath{$\displaystyle\mathrm{S}$}}_{\scriptscriptstyle
GYR}\ \mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
YP}^{\star}=\begin{bmatrix}\mathrm{m_{\scriptscriptstyle{GYR,11}}}&\mathrm{m_{\scriptscriptstyle{GYR,12}}}&\mathrm{m_{\scriptscriptstyle{GYR,13}}}\\\
\mathrm{m_{\scriptscriptstyle{GYR,21}}}&\mathrm{m_{\scriptscriptstyle{GYR,22}}}&\mathrm{m_{\scriptscriptstyle{GYR,23}}}\\\
\mathrm{m_{\scriptscriptstyle{GYR,31}}}&\mathrm{m_{\scriptscriptstyle{GYR,32}}}&\mathrm{m_{\scriptscriptstyle{GYR,33}}}\end{bmatrix}}$
(47) $\displaystyle\displaystyle\mathrm{\approx}$
$\displaystyle\displaystyle\mathrm{\begin{bmatrix}\mathrm{s_{\scriptscriptstyle{GYR,1}}}&\mathrm{\alpha_{\scriptscriptstyle{GYR,23}}\
\left(s_{\scriptscriptstyle{GYR,1}}-s_{\scriptscriptstyle{GYR,2}}\right)}&\mathrm{\alpha_{\scriptscriptstyle{GYR,32}}\
\left(s_{\scriptscriptstyle{GYR,3}}-s_{\scriptscriptstyle{GYR,1}}\right)}\\\
\mathrm{\alpha_{\scriptscriptstyle{GYR,13}}\
\left(s_{\scriptscriptstyle{GYR,1}}-s_{\scriptscriptstyle{GYR,2}}\right)}&\mathrm{s_{\scriptscriptstyle{GYR,2}}}&\mathrm{\alpha_{\scriptscriptstyle{GYR,31}}\
\left(s_{\scriptscriptstyle{GYR,2}}-s_{\scriptscriptstyle{GYR,3}}\right)}\\\
\mathrm{\alpha_{\scriptscriptstyle{GYR,12}}\
\left(s_{\scriptscriptstyle{GYR,3}}-s_{\scriptscriptstyle{GYR,1}}\right)}&\mathrm{\alpha_{\scriptscriptstyle{GYR,21}}\
\left(s_{\scriptscriptstyle{GYR,2}}-s_{\scriptscriptstyle{GYR,3}}\right)}&\mathrm{s_{\scriptscriptstyle{GYR,3}}}\end{bmatrix}}$
$\displaystyle\displaystyle\mathrm{m_{\scriptscriptstyle{GYR,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(1,\
s_{\scriptscriptstyle GYR}^{2}\right)\hskip 122.0pti=j}$ (48)
$\displaystyle\displaystyle\mathrm{m_{\scriptscriptstyle{GYR,ij}}}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\mathrm{N\left(0,\ 2\
\alpha_{\scriptscriptstyle GYR}^{2}\ s_{\scriptscriptstyle
GYR}^{2}\right)=N\left(0,\ m_{\scriptscriptstyle GYR}^{2}\right)\hskip
20.0pti\neq j}$ (49)

###  3.8 Mounting of Inertial Sensors

Equations (43) and (46) contain the relationships between the specific force
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle P}}$ and inertial angular velocity
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle P}}$ and their measurements
$\displaystyle\mathrm{\left(\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle P},\
\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IP}^{\scriptscriptstyle P}}\right)}$, when evaluated and viewed in the
platform frame $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$. However, from
the point of view of the navigation system, both magnitudes need to be
evaluated and viewed in the body frame
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ instead of
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$. These equations thus need to
be modified so they relate
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$ with
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$ as well as
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$ with
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$, respectively, as described in section 3.9 below.
To do that, it is necessary to define the relative pose (position plus
attitude) between the $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ and
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ frames. Note that the IMU,
represented by the platform frame $\displaystyle\mathrm{F_{\scriptscriptstyle
P}}$, should be mounted in the aircraft as close as possible to the center of
gravity (this reduces errors, as described in section 3.9), and loosely
aligned with the aircraft body axes.

To provide realism to the stochastic simulation, this article assumes that the
real displacement
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}$ between the two frames is deterministic, while
the relative rotation
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}^{\scriptscriptstyle
BP}=\left[\psi^{\scriptscriptstyle P},\ \theta^{\scriptscriptstyle P},\
\xi^{\scriptscriptstyle P}\right]^{T}}$ is stochastic. In this way, each
simulated aircraft represented by its specific aircraft seed
$\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$ exhibits a slightly
different IMU platform attitude with respect to the aircraft body:

  * •

As the IMU reference point is fixed with respect to the structure but the
aircraft center of mass slowly moves as the fuel load diminishes, it is
possible to establish a linear model that provides the displacement between
the origins of both frames according to the aircraft mass202020The aircraft
masses $\displaystyle\mathrm{m_{full}}$ and $\displaystyle\mathrm{m_{empty}}$
when the fuel tank is fully loaded or empty are inputs, as are the
displacements between the IMU reference point and the aircraft center of mass
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{{\scriptscriptstyle
BP},full}^{\scriptscriptstyle B}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{{\scriptscriptstyle
BP},empty}^{\scriptscriptstyle B}}$.:

$\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}=\mbox{\boldmath{$\displaystyle\mathrm{f}$}}\left(m\right)=\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{{\scriptscriptstyle
BP},full}^{\scriptscriptstyle
B}+\dfrac{m_{full}-m}{m_{full}-m_{empty}}\;\left(\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{{\scriptscriptstyle
BP},empty}^{\scriptscriptstyle
B}-\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{{\scriptscriptstyle
BP},full}^{\scriptscriptstyle B}\right)}$ (50)

  * •

The platform Euler angles respond to the stochastic model provided by (51), in
which each specific Euler angle is obtained as the product of the standard
deviations ($\displaystyle\mathrm{\mathrm{\sigma_{\psi^{\scriptscriptstyle
P}}}}$, $\displaystyle\mathrm{\mathrm{\sigma_{\theta^{\scriptscriptstyle
P}}}}$, $\displaystyle\mathrm{\mathrm{\sigma_{\xi^{\scriptscriptstyle P}}}}$)
by a single realization of a standard normal random variable
$\displaystyle\mathrm{N\left(0,\,1\right)}$
($\displaystyle\mathrm{N_{\psi^{\scriptscriptstyle P}}}$,
$\displaystyle\mathrm{N_{\theta^{\scriptscriptstyle P}}}$, and
$\displaystyle\mathrm{N_{\xi^{\scriptscriptstyle P}}}$).

$\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}^{\scriptscriptstyle
BP}=\left[\mathrm{\sigma_{\psi^{\scriptscriptstyle
P}}}\,N_{\psi^{\scriptscriptstyle P}},\
\mathrm{\sigma_{\theta^{\scriptscriptstyle P}}}\,N_{\theta^{\scriptscriptstyle
P}},\ \mathrm{\sigma_{\xi^{\scriptscriptstyle P}}}N_{\xi^{\scriptscriptstyle
P}}\right]^{T}}$ (51)

Concept | Variable | Value | Unit
---|---|---|---
True Yaw Error | $\displaystyle\mathrm{\mathrm{\sigma_{\psi^{\scriptscriptstyle P}}}}$ | $\displaystyle\mathrm{0.5}$ | [$\displaystyle\mathrm{deg}$]
True Pitch Error | $\displaystyle\mathrm{\mathrm{\sigma_{\theta^{\scriptscriptstyle P}}}}$ | $\displaystyle\mathrm{2.0}$ | [$\displaystyle\mathrm{deg}$]
True Bank Error | $\displaystyle\mathrm{\mathrm{\sigma_{\xi^{\scriptscriptstyle P}}}}$ | $\displaystyle\mathrm{0.1}$ | [$\displaystyle\mathrm{deg}$]
Platform Position Estimation Error | $\displaystyle\mathrm{\mathrm{\sigma_{\hat{T}_{\scriptscriptstyle BP}^{\scriptscriptstyle B}}}}$ | $\displaystyle\mathrm{0.01}$ | [$\displaystyle\mathrm{m}$]
Platform Angular Estimation Error | $\displaystyle\mathrm{\mathrm{\sigma_{\hat{\phi}^{\scriptscriptstyle BP}}}}$ | $\displaystyle\mathrm{0.03}$ | [$\displaystyle\mathrm{deg}$]
Table 7: IMU mounting accuracy values

Once the real pose between the $\displaystyle\mathrm{F_{\scriptscriptstyle
P}}$ and $\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ frames is
established, it is necessary to specify its estimation employed by the IMU in
the comprehensive model introduced in section 3.9, which is discussed in
section 6.2. Stochastic models are employed in the simulation for both the
translation
$\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}$ and rotation
$\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}}^{\scriptscriptstyle
BP}}$, changing their values from one aircraft seed
$\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$ to the next:

$\displaystyle\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}+\left[\mathrm{\sigma_{\hat{T}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}}\,N_{\hat{T}_{\scriptscriptstyle
BP,1}^{\scriptscriptstyle P}},\ \mathrm{\sigma_{\hat{T}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}}\,N_{\hat{T}_{\scriptscriptstyle
BP,2}^{\scriptscriptstyle P}},\ \mathrm{\sigma_{\hat{T}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}}N_{\hat{T}_{\scriptscriptstyle
BP,3}^{\scriptscriptstyle P}}\right]^{T}}$ (52)
$\displaystyle\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}}^{\scriptscriptstyle
BP}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}^{\scriptscriptstyle
BP}+\left[\mathrm{\sigma_{\hat{\phi}^{\scriptscriptstyle
BP}}}\,N_{\hat{\psi}^{\scriptscriptstyle P}},\
\mathrm{\sigma_{\hat{\phi}^{\scriptscriptstyle
BP}}}\,N_{\hat{\theta}^{\scriptscriptstyle P}},\
\mathrm{\sigma_{\hat{\phi}^{\scriptscriptstyle
BP}}}N_{\hat{\xi}^{\scriptscriptstyle P}}\right]^{T}}$ (53)

where the default standard deviations
$\displaystyle\mathrm{\mathrm{\sigma_{\hat{T}_{\scriptscriptstyle
BP}^{\scriptscriptstyle B}}}}$ and
$\displaystyle\mathrm{\mathrm{\sigma_{\hat{\phi}^{\scriptscriptstyle BP}}}}$
are shown in table 3.8, and
$\displaystyle\mathrm{N_{\hat{\psi}^{\scriptscriptstyle P}}}$,
$\displaystyle\mathrm{N_{\hat{\theta}^{\scriptscriptstyle P}}}$,
$\displaystyle\mathrm{N_{\hat{\xi}^{\scriptscriptstyle P}}}$,
$\displaystyle\mathrm{N_{\hat{T}_{\scriptscriptstyle BP,1}^{\scriptscriptstyle
P}}}$, $\displaystyle\mathrm{N_{\hat{T}_{\scriptscriptstyle
BP,2}^{\scriptscriptstyle P}}}$,
$\displaystyle\mathrm{N_{\hat{T}_{\scriptscriptstyle BP,3}^{\scriptscriptstyle
P}}}$ are six realizations of a standard normal random variable
$\displaystyle\mathrm{N\left(0,\,1\right)}$. Note that table 3.8 lists the
default values employed in the simulation, which can be adjusted by the user.

$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}}$ can be considered quasi stationary as it slowly varies based on the
aircraft mass, and the relative position of their axes
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\phi}$}}^{\scriptscriptstyle
BP}}$ remains constant because the IMU is rigidly attached to the aircraft
structure. Although Euler angles have been employed in this section, from this
point on it is more practical to employ the rotation matrix
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
BP}}$ to represent the rotation between two different frames [21]. The time
derivatives of
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
BP}}$ are hence zero:

$\mathrm{\dot{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}}_{\scriptscriptstyle
BP}=\dot{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle
BP}=\mbox{\boldmath{$\displaystyle\mathrm{0}$}}\ \longrightarrow\
\mbox{\boldmath{$\displaystyle\mathrm{v}$}}_{\scriptscriptstyle
BP}=\mbox{\boldmath{$\displaystyle\mathrm{a}$}}_{\scriptscriptstyle
BP}=\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
BP}=\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}_{\scriptscriptstyle
BP}=\mbox{\boldmath{$\displaystyle\mathrm{0}$}}}$ (54)

###  3.9 Comprehensive Inertial Sensor Error Model

Two considerations are required to establish the measurement equations for the
inertial sensors viewed in the body frame
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$. First, let’s apply the
composition rules of section A.6 considering
$\displaystyle\mathrm{F_{\scriptscriptstyle I}}$ as
$\displaystyle\mathrm{F_{\scriptscriptstyle 0}}$,
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ as
$\displaystyle\mathrm{F_{\scriptscriptstyle 1}}$, and
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ as
$\displaystyle\mathrm{F_{\scriptscriptstyle 2}}$, which results in:

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IP}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IB}}$ (55)
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}_{\scriptscriptstyle
IP}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}_{\scriptscriptstyle
IB}}$ (56)
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}_{\scriptscriptstyle
IP}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{v}$}}_{\scriptscriptstyle
IB}+\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}\ \mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle BP}}$
(57)
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{a}$}}_{\scriptscriptstyle
IP}}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{a}$}}_{\scriptscriptstyle
IB}+\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}}_{\scriptscriptstyle
IB}\ \mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}+\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}\
\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}\ \mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle BP}}$
(58)

Second, it is also necessary to consider that as
$\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle
BP}}$ is a rotation matrix in which all rows and columns are unitary vectors,
the projection of the $\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ frame
bias and white noise errors
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle P}}$ and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle P}}$ onto the
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ frame does not change their
stochastic properties:

$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle B}\left(s\,\Delta t\right)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle
BP}\ \mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle
P}=B_{0\scriptscriptstyle{ACC}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{u0,\scriptscriptstyle
ACC}+\sigma_{u\scriptscriptstyle{ACC}}\,\Delta
t^{1/2}\,\sum_{i=1}^{s}\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{ui,\scriptscriptstyle
ACC}+\dfrac{\sigma_{v\scriptscriptstyle{ACC}}}{\Delta
t^{1/2}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{vs,\scriptscriptstyle
ACC}}$ (59)
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle B}\left(s\,\Delta t\right)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle
BP}\ \mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle
P}=B_{0\scriptscriptstyle{GYR}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{u0,\scriptscriptstyle
GYR}+\sigma_{u\scriptscriptstyle{GYR}}\,\Delta
t^{1/2}\,\sum_{i=1}^{s}\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{ui,\scriptscriptstyle
GYR}+\dfrac{\sigma_{v\scriptscriptstyle{GYR}}}{\Delta
t^{1/2}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{vs,\scriptscriptstyle
GYR}}$ (60)

As the inertial angular velocity does not change when evaluated in the
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ and
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ frames (55), its measurement
in the body frame can be derived from (46) by first projecting it from
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ to
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ based on the real rotation
matrix
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
BP}}$ and then projecting back the measurement into
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ based on the estimated
rotation matrix
$\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle
BP}}$. The bias and white noise error is also projected according to (60):

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B}=\hat{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle BP}\
\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle GYR}\
{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle BP}}^{T}\
\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,GYR}^{\scriptscriptstyle B}}$ (61)

The expression for the specific force measurement is significantly more
complex because the back and forth transformations of the specific force
between the $\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ and
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ frames need to consider the
influence of the lever arm
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}}$, as indicated in (58). The additional terms introduce errors in the
measurements, so as indicated in section 3.8 it is desirable to locate the IMU
as close as possible to the aircraft center of mass.

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B}=\hat{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle BP}\
\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle ACC}\
{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle BP}}^{T}\
\left(\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B}+\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}+\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}\right)-\widehat{\hat{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\hat{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}-\widehat{\hat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\widehat{\hat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\hat{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle B}}$ (62)

Note that this expression can not be directly evaluated as the estimated
values for the inertial angular velocity and acceleration
($\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B},\
\hat{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$) are unknown by the IMU until obtained by the
navigation filter. The IMU can however rely on the gyroscope readings,
directly replacing
$\displaystyle\mathrm{\hat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$ with
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$ and computing
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$ based on the difference between the present and
previous
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}}$ readings, resulting in:

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B}=\hat{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}}_{\scriptscriptstyle BP}\
\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle ACC}\
{\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle BP}}^{T}\
\left(\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle
B}+\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}+\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\widehat{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}\right)-\widehat{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\alpha}$}}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\hat{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}-\widehat{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\widehat{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\
\hat{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}}_{\scriptscriptstyle
BP}^{\scriptscriptstyle
B}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
BW,ACC}^{\scriptscriptstyle B}}$ (63)

Table 3.9 lists the error sources contained in the comprehensive inertial
sensor error model represented by (61, 63). The first two columns list the
different error sources, while the third column specifies their origin
according to the criterion established in the first paragraph of section 3.1.
The section where each error is described appears on the fourth column,
followed by the seeds employed to ensure the results variability for different
aircraft ($\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$) as well
as different flights ($\displaystyle\mathrm{\Upsilon_{j,\scriptscriptstyle
F}}$):

Error | Main Source | Description | Seeds
---|---|---|---
Bias Offset | $\displaystyle\mathrm{B_{0\scriptscriptstyle{ACC}},\ B_{0\scriptscriptstyle{GYR}}}$ | run-to-run | 3.2 |  $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle ACC}}$, $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle GYR}}$ | $\displaystyle\mathrm{\Upsilon_{j,\scriptscriptstyle F}}$
Bias Drift | $\displaystyle\mathrm{\sigma_{u\scriptscriptstyle{ACC}},\ \sigma_{u\scriptscriptstyle{GYR}}}$ | in-run | 3.2 |  $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle ACC}}$, $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle GYR}}$ | $\displaystyle\mathrm{\Upsilon_{j,\scriptscriptstyle F}}$
System Noise | $\displaystyle\mathrm{\sigma_{v\scriptscriptstyle{ACC}},\ \sigma_{v\scriptscriptstyle{GYR}}}$ | in-run | 3.2 |  $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle ACC}}$, $\displaystyle\mathrm{\upsilon_{j,\scriptscriptstyle F,\scriptscriptstyle GYR}}$ | $\displaystyle\mathrm{\Upsilon_{j,\scriptscriptstyle F}}$
Scale Factor | $\displaystyle\mathrm{s_{\scriptscriptstyle ACC},\ s_{\scriptscriptstyle GYR}}$ | fixed & T |  3.6, 3.7 |  $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle ACC}}$, $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle GYR}}$ | $\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$
Cross Coupling | $\displaystyle\mathrm{m_{\scriptscriptstyle ACC},\ m_{\scriptscriptstyle GYR}}$ | fixed |  3.6, 3.7 |  $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle ACC}}$, $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle GYR}}$ | $\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$
Lever Arm | $\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle BP},\ \mathrm{\sigma_{\hat{T}_{\scriptscriptstyle BP}^{\scriptscriptstyle B}}}}$ | fixed | 3.8 | $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle PLAT}}$ | $\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$
IMU Attitude | $\displaystyle\mathrm{\mathrm{\sigma_{\psi^{\scriptscriptstyle P}}},\ \mathrm{\sigma_{\theta^{\scriptscriptstyle P}}},\ \mathrm{\sigma_{\xi^{\scriptscriptstyle P}}},\ \mathrm{\sigma_{\hat{\phi}^{\scriptscriptstyle BP}}}}$ | fixed | 3.8 | $\displaystyle\mathrm{\upsilon_{i,\scriptscriptstyle A,\scriptscriptstyle PLAT}}$ | $\displaystyle\mathrm{\Upsilon_{i,\scriptscriptstyle A}}$
Table 8: Inertial sensor error sources

Note that $\displaystyle\mathrm{B_{0\scriptscriptstyle{ACC}}}$,
$\displaystyle\mathrm{B_{0\scriptscriptstyle{GYR}}}$,
$\displaystyle\mathrm{\sigma_{u\scriptscriptstyle{ACC}}}$,
$\displaystyle\mathrm{\sigma_{u\scriptscriptstyle{GYR}}}$,
$\displaystyle\mathrm{\sigma_{v\scriptscriptstyle{ACC}}}$,
$\displaystyle\mathrm{\sigma_{v\scriptscriptstyle{GYR}}}$,
$\displaystyle\mathrm{s_{\scriptscriptstyle ACC}}$,
$\displaystyle\mathrm{s_{\scriptscriptstyle GYR}}$,
$\displaystyle\mathrm{m_{\scriptscriptstyle ACC}}$, and
$\displaystyle\mathrm{m_{\scriptscriptstyle GYR}}$ are taken from the right
hand column of tables 3.10 and 3.10 (section 3.10), while the remaining error
sources come from section 3.8. It is worth pointing out that all errors are
modeled as stochastic variables or processes (with the exception of the
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{T}$}}_{\scriptscriptstyle
BP}}$ displacement between the body and platform frames, which is
deterministic), as expressions (61, 63) rely on the errors provided by (59,
60), the scale and cross coupling matrices given by (34, 47), and the
transformations given by (50, 51, 52, 53).

In the case of the accelerometer triad, the stochastic nature of the fixed and
run-to-run error contributions to the model relies on three realizations of
normal distributions for the bias offset, three for the scale factor errors,
three for the cross coupling errors, and nine for the mounting errors, while
the in-run error contributions require three realizations each for the bias
drift and system noise at every discrete sensor measurement. The gyroscope
triad is similar but requires six realizations to model the cross coupling
errors instead of three, while using the same six realizations that the
accelerometer triad to model the true and estimated rotation between the
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$ and
$\displaystyle\mathrm{F_{\scriptscriptstyle P}}$ frames.

Expressions (61, 63) can be rewritten to show the measurements as functions of
the full errors
($\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{E}$}}_{\scriptscriptstyle
ACC},\,\mbox{\boldmath{$\displaystyle\mathrm{E}$}}_{\scriptscriptstyle
GYR}}$), which represent all the errors introduced by the inertial sensors
with the exception of white noise.

$\displaystyle\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\left(s\,\Delta t\right)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{f}$}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\left(s\,\Delta
t\right)+\mbox{\boldmath{$\displaystyle\mathrm{E}$}}_{\scriptscriptstyle
ACC}\left(s\,\Delta t\right)+\dfrac{\sigma_{v\scriptscriptstyle{ACC}}}{\Delta
t^{1/2}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{vs,\scriptscriptstyle
ACC}}$ (64)
$\displaystyle\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\left(s\,\Delta t\right)}$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{\omega}$}}_{\scriptscriptstyle
IB}^{\scriptscriptstyle B}\left(s\,\Delta
t\right)+\mbox{\boldmath{$\displaystyle\mathrm{E}$}}_{\scriptscriptstyle
GYR}\left(s\,\Delta t\right)+\dfrac{\sigma_{v\scriptscriptstyle{GYR}}}{\Delta
t^{1/2}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{vs,\scriptscriptstyle
GYR}}$ (65)

###  3.10 Baseline Inertial Sensors

The value of all parameters present in the previous models can be specified by
the user when employing the simulation [1, 2]. This article however provides
default values that represent the inertial sensors that could nowadays be
installed onboard a small size fixed wing low SWaP autonomous aircraft, and
which are referred to as the baseline configuration. In the case of the
gyroscopes, they correspond to the MEMS gyroscopes installed inside the Analog
Devices ADIS16488A IMU [22]. Table 3.10 shows its performances, which have
been taken from the data sheet when possible, and corrected when suspicious. A
calibration process as that described in section 6.1 is assumed to eliminate
$\displaystyle\mathrm{95\%}$ of the scale factor and cross coupling errors.

GYR Baseline | Spec | Unit | Variable | Value | Calibration | Unit
---|---|---|---|---|---|---
In-Run Bias Stability (1 $\displaystyle\mathrm{\sigma}$) | 5.10 | [$\displaystyle\mathrm{deg/hr}$] | $\displaystyle\mathrm{\sigma_{u\scriptscriptstyle{GYR}}}$ | $\displaystyle\mathrm{1.42\cdot 10^{-4}}$ | $\displaystyle\mathrm{1.42\cdot 10^{-4}}$ | [$\displaystyle\mathrm{deg/sec^{1.5}}$]
Angle Random Walk (1 $\displaystyle\mathrm{\sigma}$) | 0.26 | [$\displaystyle\mathrm{deg/hr^{0.5}}$] | $\displaystyle\mathrm{\sigma_{v\scriptscriptstyle{GYR}}}$ | $\displaystyle\mathrm{4.30\cdot 10^{-3}}$ | $\displaystyle\mathrm{4.30\cdot 10^{-3}}$ | [$\displaystyle\mathrm{deg/sec^{0.5}}$]
Nonlinearity212121The $\displaystyle\mathrm{0.01\%}$ scale factor error obtained in [22] is considered too optimistic and hence modified to $\displaystyle\mathrm{0.03\%=3.00\cdot 10^{-4}\ [-]}$. | 0.01 | [%] | $\displaystyle\mathrm{s_{\scriptscriptstyle GYR}}$ | $\displaystyle\mathrm{3.00\cdot 10^{-4}}$ | $\displaystyle\mathrm{1.50\cdot 10^{-5}}$ | [-]
Misalignment |  $\displaystyle\mathrm{\pm}$ 0.05 | [$\displaystyle\mathrm{deg}$] | $\displaystyle\mathrm{m_{\scriptscriptstyle GYR}}$ | $\displaystyle\mathrm{8.70\cdot 10^{-4}}$ | $\displaystyle\mathrm{4.35\cdot 10^{-5}}$ | [-]
Bias Repeatability (1 $\displaystyle\mathrm{\sigma}$) |  $\displaystyle\mathrm{\pm}$ 0.2 | [$\displaystyle\mathrm{deg/sec}$] | $\displaystyle\mathrm{B_{0\scriptscriptstyle{GYR}}}$ | $\displaystyle\mathrm{2.00\cdot 10^{-1}}$ | $\displaystyle\mathrm{2.00\cdot 10^{-1}}$ | [$\displaystyle\mathrm{deg/sec}$]
Table 9: Performance of “Baseline” gyroscopes

The previous table show the performances in three columns of data. The left
most column ("Spec") corresponds to data taken directly from the sensors
specifications, which get converted in the middle column ("Value") to the
parameters shown in sections 3.2 through 3.9222222The conversion between bias
instability and $\displaystyle\mathrm{\sigma_{u}}$ uses a period of 100
[$\displaystyle\mathrm{sec}$], as noted in section 3.4.. The right column
("Calibration") contains the specifications employed in the simulation after
the calibration process, which reduces the scale factor and cross coupling
errors by $\displaystyle\mathrm{95\%}$.

ACC Baseline | Spec | Unit | Variable | Value | Calibration | Unit
---|---|---|---|---|---|---
In-Run Bias Stability (1 $\displaystyle\mathrm{\sigma}$) | 0.07 | [$\displaystyle\mathrm{mg}$] | $\displaystyle\mathrm{\sigma_{u\scriptscriptstyle{ACC}}}$ | $\displaystyle\mathrm{6.86\cdot 10^{-5}}$ | $\displaystyle\mathrm{6.86\cdot 10^{-5}}$ | [$\displaystyle\mathrm{m/sec^{2.5}}$]
Velocity Random Walk (1 $\displaystyle\mathrm{\sigma}$) | 0.029 | [$\displaystyle\mathrm{m/sec/hr^{0.5}}$] | $\displaystyle\mathrm{\sigma_{v\scriptscriptstyle{ACC}}}$ | $\displaystyle\mathrm{4.83\cdot 10^{-4}}$ | $\displaystyle\mathrm{4.83\cdot 10^{-4}}$ | [$\displaystyle\mathrm{m/sec^{1.5}}$]
Nonlinearity | 0.1 | [%] | $\displaystyle\mathrm{s_{\scriptscriptstyle ACC}}$ | $\displaystyle\mathrm{1.00\cdot 10^{-3}}$ | $\displaystyle\mathrm{5.00\cdot 10^{-5}}$ | [-]
Misalignment |  $\displaystyle\mathrm{\pm}$ 0.035 | [$\displaystyle\mathrm{deg}$] | $\displaystyle\mathrm{m_{\scriptscriptstyle ACC}}$ | $\displaystyle\mathrm{6.11\cdot 10^{-4}}$ | $\displaystyle\mathrm{3.05\cdot 10^{-5}}$ | [-]
Bias Repeatability (1 $\displaystyle\mathrm{\sigma}$) |  $\displaystyle\mathrm{\pm}$ 16 | [$\displaystyle\mathrm{mg}$] | $\displaystyle\mathrm{B_{0\scriptscriptstyle{ACC}}}$ | $\displaystyle\mathrm{1.57\cdot 10^{-1}}$ | $\displaystyle\mathrm{1.57\cdot 10^{-1}}$ | [$\displaystyle\mathrm{m/sec^{2}}$]
Table 10: Performance of “Baseline” accelerometers

The baseline accelerometers also correspond to the MEMS accelerometers
installed inside the Analog Devices ADIS16488A IMU [22]. All the
specifications shown in table 3.10 have been taken from the data sheet. As in
the case of the gyroscopes, a calibration process as that described in section
6.1 is assumed to eliminate $\displaystyle\mathrm{95\%}$ of the scale factor
and cross coupling errors232323The conversion between bias instability and
$\displaystyle\mathrm{\sigma_{u}}$ uses a period of 100
[$\displaystyle\mathrm{sec}$], as noted in section 3.4..

##  4 Non Inertial Sensors

This section describes the different non-inertial sensors usually installed
onboard a fixed wind autonomous aircraft, such as a triad of _magnetometers_
to measure the Earth magnetic field, a _GNSS receiver_ that provides absolute
position and velocity measurements, and the _air data system_ , which in
addition to the pressure altitude and temperature also provides a measurement
of the airspeed and the airflow angles.

###  4.1 Magnetometers

Magnetometers measure magnetic field intensity along a given direction and are
very useful for estimating the aircraft heading. Although other types exist,
magnetoinductive and magnetoresistive sensors are generally employed for
navigation due to their accuracy and small size [5, 7]. As with the inertial
sensors, three orthogonal magnetometers are usually employed in a strapdown
configuration to measure the magnetic field with respect to the body frame
$\displaystyle\mathrm{F_{\scriptscriptstyle B}}$.

Unfortunately magnetometers do not only measure the Earth magnetic field
$\displaystyle\mathrm{B}$, but also that generated by the aircraft permanent
magnets and electrical equipment (known as hard iron magnetism), as well as
the magnetic field disturbances generated by the aircraft ferrous materials
(soft iron magnetism). For that reason, the magnetometers should be placed in
a location inside the aircraft that minimizes these errors. On the positive
side, magnetometers do not exhibit the bias instability present in inertial
sensors, and the error of an individual sensor can be properly modeled by the
combination of bias offset and white noise. A triad of magnetometers capable
of measuring the magnetic field in three directions adds the same scale factor
(nonlinearity) and cross coupling (misalignment) errors as those present in
the inertial sensors, together with the transformation between the magnetic
axes and the body ones.

Modeling the behavior of a triad of magnetometers is simpler but less precise
than that of inertial sensors, as the effect of the fixed hard iron magnetism
is indistinguishable from that of the run-to-run bias offset, while the fixed
effect of soft iron magnetism is indistinguishable from that of the scale
factor and cross coupling error matrix. This has several consequences. First
of all is that magnetometers can not be calibrated at the laboratory before
being mounted in the aircraft as in the case of inertial sensors (section
6.1), but are instead calibrated once attached to the aircraft by a process
known as swinging (section 6.3), which is less precise as the aircraft
attitude during swinging can not be determined with so much accuracy as it
would be in a laboratory setting. Second is that defining a magnetic platform
frame to then transform the results into body axes serves no purpose, as the
magnetometer readings are only valid, this is, contain the effects of hard and
soft iron magnetism, once they are attached to the aircraft, and then they can
be directly measured in body axes. And third is that percentage wise the
errors induced by the magnetometers are bigger than those of the inertial
sensors. The simulation relies on the following model:

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}}^{\scriptscriptstyle
B}=\mbox{\boldmath{$\displaystyle\mathrm{B}$}}_{0,\scriptscriptstyle
MAG}+\mbox{\boldmath{$\displaystyle\mathrm{B}$}}_{\scriptscriptstyle{HI,MAG}}+\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
MAG}\,\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
BN}\,\mbox{\boldmath{$\displaystyle\mathrm{B}$}}^{\scriptscriptstyle
N}+\mbox{\boldmath{$\displaystyle\mathrm{e}$}}_{\scriptscriptstyle
W,MAG}^{\scriptscriptstyle B}}$ (66)

where
$\displaystyle\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}}^{\scriptscriptstyle
B}}$ is the measurement viewed in $\displaystyle\mathrm{F_{\scriptscriptstyle
B}}$,
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}_{0,\scriptscriptstyle
MAG}}$ is the run-to-run bias offset,
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}_{\scriptscriptstyle{HI,MAG}}}$
is the fixed hard iron magnetism,
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
MAG}}$ is a fixed matrix combining the effects of soft iron magnetism with the
scale factor and cross coupling errors, and
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}^{\scriptscriptstyle
N}}$ is the real magnetic field including local anomalies.

MAG Baseline | Spec | Unit | Variable | Value | Comp. | Swinging | Unit
---|---|---|---|---|---|---|---
Output Noise | 5 | [$\displaystyle\mathrm{nT\cdot sec^{0.5}}$] | $\displaystyle\mathrm{\sigma_{v,\scriptscriptstyle MAG}}$ | $\displaystyle\mathrm{5.00\cdot 10^{0}}$ | $\displaystyle\mathrm{5.00\cdot 10^{0}}$ | $\displaystyle\mathrm{5.00\cdot 10^{0}}$ | [$\displaystyle\mathrm{nT\cdot sec^{0.5}}$]
Nonlinearity | 0.5 | [%] | $\displaystyle\mathrm{s_{\scriptscriptstyle MAG}}$ | $\displaystyle\mathrm{5.00\cdot 10^{-3}}$ | $\displaystyle\mathrm{7.50\cdot 10^{-3}}$ | $\displaystyle\mathrm{7.50\cdot 10^{-4}}$ | [-]
Misalignment |  $\displaystyle\mathrm{\pm}$ 0.35 | [$\displaystyle\mathrm{deg}$] | $\displaystyle\mathrm{m_{\scriptscriptstyle MAG}}$ | $\displaystyle\mathrm{6.11\cdot 10^{-3}}$ | $\displaystyle\mathrm{9.16\cdot 10^{-3}}$ | $\displaystyle\mathrm{9.16\cdot 10^{-4}}$ | [-]
Bias (1 $\displaystyle\mathrm{\sigma}$) | $\displaystyle\mathrm{\pm}$ 1500 | [$\displaystyle\mathrm{nT}$] | $\displaystyle\mathrm{B_{\scriptscriptstyle{HI,MAG}}}$ | $\displaystyle\mathrm{1.50\cdot 10^{3}}$ | $\displaystyle\mathrm{1.75\cdot 10^{3}}$ | $\displaystyle\mathrm{1.75\cdot 10^{2}}$ | [$\displaystyle\mathrm{nT}$]
Repeatability | $\displaystyle\mathrm{B_{0,\scriptscriptstyle MAG}}$ | $\displaystyle\mathrm{5.00\cdot 10^{2}}$ | $\displaystyle\mathrm{5.00\cdot 10^{2}}$ | [$\displaystyle\mathrm{nT}$]
Table 11: Performance of “Baseline” magnetometers

The baseline magnetometer features employed in the simulation are shown in
table 4.1 in the same format as those of the inertial sensors, where the white
noise has been taken from [5] and the rest of the parameters correspond to the
magnetometers present in the Analog Devices ADIS16488A IMU [22]. Although the
value of hard and soft iron magnetism in aircraft is rather small, the author
has not been able to obtain trusted values for them. To avoid eliminating
sources of error, the author has decided to increase by 50% the values for
bias offset, scale factor, and cross coupling errors found in the literature,
as shown in the column named "Compensation". As both result in a similar
effect, the author expects that the realism of the results will not be
adversely affected. In the case of the bias, the author has assigned most of
the error to the fixed hard iron error
$\displaystyle\mathrm{B_{\scriptscriptstyle{HI,MAG}}}$ and the remaining to
the run-to-run bias offset $\displaystyle\mathrm{B_{0,\scriptscriptstyle
MAG}}$.

A swinging process as that described in section 6.3 is assumed to eliminate
90% of the fixed error contributions, this is, the hard iron magnetism, the
scale factor, and the cross coupling error242424The soft iron effect is
combined with the scale factor and cross coupling errors.. This number is
inferior to the 95% reduction achieved by the calibration of inertial sensors
because although the lack of magnetometer bias drift facilitates calibration,
the determination of the body attitude during calibration is inherently less
accurate working with the complete aircraft than just with the IMU in a
laboratory setting. The final results can be found in the rightmost column
above.

The final model for a triaxial magnetometer thus includes contributions from
the hard iron magnetism, bias offset, system noise, soft iron magnetism, scale
factor, and cross coupling errors:

$\mathrm{\widetilde{\mbox{\boldmath{$\displaystyle\mathrm{B}$}}}^{\scriptscriptstyle
B}\left(s\,\Delta
t\right)=B_{\scriptscriptstyle{HI,MAG}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{\scriptscriptstyle{HI},\scriptscriptstyle
MAG}+B_{0,\scriptscriptstyle
MAG}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{u0,\scriptscriptstyle
MAG}+\mbox{\boldmath{$\displaystyle\mathrm{M}$}}_{\scriptscriptstyle
MAG}\,\mbox{\boldmath{$\displaystyle\mathrm{R}$}}_{\scriptscriptstyle
BN}\,\mbox{\boldmath{$\displaystyle\mathrm{B}$}}_{\scriptscriptstyle
REAL}^{\scriptscriptstyle N}+\dfrac{\sigma_{v,\scriptscriptstyle MAG}}{\Delta
t^{1/2}}\,\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{vs,\scriptscriptstyle
MAG}}$ (67)

where
$\displaystyle\mathrm{\mbox{\boldmath{$\displaystyle\mathrm{N}$}}_{\scriptscriptstyle{HI},\scriptscriptstyle