Using Inertial Sensors for Position and Orientation Estimation

TL;DR

This paper reviews methods for estimating position and orientation using MEMS inertial sensors, highlighting algorithms like Kalman filters and smoothing techniques, with evaluations on experimental and simulated data.

Contribution

It provides a comprehensive overview of signal processing algorithms for inertial sensor-based estimation, emphasizing modeling choices and practical implementations.

Findings

Algorithms like extended Kalman filter and complementary filter effectively estimate position and orientation.

Signal processing techniques can mitigate drift in inertial sensor measurements.

Experimental and simulated data validate the effectiveness of discussed algorithms.

Abstract

In recent years, MEMS inertial sensors (3D accelerometers and 3D gyroscopes) have become widely available due to their small size and low cost. Inertial sensor measurements are obtained at high sampling rates and can be integrated to obtain position and orientation information. These estimates are accurate on a short time scale, but suffer from integration drift over longer time scales. To overcome this issue, inertial sensors are typically combined with additional sensors and models. In this tutorial we focus on the signal processing aspects of position and orientation estimation using inertial sensors. We discuss different modeling choices and a selected number of important algorithms. The algorithms include optimization-based smoothing and filtering as well as computationally cheaper extended Kalman filter and complementary filter implementations. The quality of their estimates is illustrated using both experimental and simulated data.

Tables

Table 1. Table 4.1 : Mean RMSE values over 100 100 100 Monte Carlo simulations estimating orientation using only inertial measurements.

RMSE		Roll [∘]	Pitch [∘]	Heading [∘]
Smoothing optimization (Alg. 1)	0.39	0.39	7.46
Filtering optimization (Alg. 2)	0.46	0.46	7.46
EKF quaternions (Alg. 3)	0.46	0.46	17.49
EKF orientation deviation (Alg. 4)	0.46	0.46	7.46

Table 2. Table 4.2: RMSE of the orientation estimates obtained using Algorithms  1 – 5 and the experimental data presented in Figure  4.3 .

RMSE		Roll [∘]	Pitch [∘]	Heading [∘]
Smoothing optimization (Alg. 1)	1.03	0.48	0.81
Filtering optimization (Alg. 2)	1.14	0.56	1.28
EKF quaternions (Alg. 3)	1.13	0.52	1.04
EKF orientation deviation (Alg. 4)	1.14	0.56	1.28
Complementary filter (Alg. 5)	0.95	0.62	1.55

Table 3. Table 4.3: Mean RMSE of the orientation estimates from 100 100 100 Monte Carlo simulations using Algorithms  1 – 5 .

RMSE		Roll [∘]	Pitch [∘]	Heading [∘]
Smoothing optimization (Alg. 1)	0.39	0.39	2.30
Filtering optimization (Alg. 2)	0.45	0.45	3.54
EKF quaternions (Alg. 3)	0.45	0.45	3.57
EKF orientation deviation (Alg. 4)	0.45	0.45	3.55
Complementary filter (Alg. 5), $\alpha =0.07$	1.44	1.43	4.39
Complementary filter (Alg. 5), $\alpha =0.7$	0.47	0.47	12.98

Table 4. Table 4.4: Mean RMSE of the orientation estimates from 100 100 100 Monte Carlo simulations. The estimate of the initial orientation is assumed to be normal distributed around the true initial orientation with a standard deviation of 20 ∘ superscript 20 20^{\circ} .

RMSE		Roll [∘]	Pitch [∘]	Heading [∘]
Smoothing optimization (Alg. 1)	0.39	0.39	2.29
Filtering optimization (Alg. 2)	1.06	0.95	3.55
EKF quaternions (Alg. 3)	1.08	0.97	4.41
EKF orientation deviation (Alg. 4)	1.08	0.96	3.57
Complementary filter (Alg. 5), $\alpha =0.07$	3.02	2.77	4.48
Complementary filter (Alg. 5), $\alpha =0.7$	1.13	1.01	12.99

Table 5. Table 4.5: Mean RMSE of the position and orientation estimates from 100 100 100 Monte Carlo simulations. Considered are a stationary sensor, a sensor with constant acceleration and two cases of time-varying accelerations with different magnitudes.

RMSE		Roll [∘]	Pitch [∘]	Heading [∘]	Position [ c m ]
Stationary	0.41	0.41	11.84	0.97
Constant acceleration	1.23	0.46	11.37	0.97
Acceleration	0.41	0.41	2.68	0.97
${\left(y_{\text{a},t}\right)}_y\sim 𝒩(0,0.5)$
Acceleration	0.46	0.39	0.87	0.97
${\left(y_{\text{a},t}\right)}_y\sim 𝒩(0,5)$

Table 6. Table 5.1: Mean RMSE of the orientation estimates from 100 100 100 Monte Carlo simulations in the presence of a gyroscope bias that is being estimated.

RMSE		Roll [∘]	Pitch [∘]	Heading [∘]
Smoothing optimization	0.39	0.39	2.29
EKF orientation deviation	0.46	0.46	4.20

Table 7. Table 5.2: Mean and standard deviation of the gyroscope estimates over 500 Monte Carlo simulations with ( 0.05   0.01 − 0.04 ) 𝖳 ​ \radianpersecond superscript 0.050.01 0.04 𝖳 \radianpersecond \left(0.05\,\,\,0.01\,\,-0.04\right)^{\mathsf{T}}~{}\radianpersecond . Considered are the cases of ML estimation and MAP estimation by including the gyroscope bias as an unknown in a smoothing algorithm with a prior on the gyroscope bias of δ ω ∼ 𝒩 ​ ( 0 , σ δ ω 2 ​ ℐ 3 ) similar-to subscript 𝛿 𝜔 𝒩 0 subscript superscript 𝜎 2 subscript 𝛿 𝜔 subscript ℐ 3 \delta_{\omega}\sim\mathcal{N}(0,\sigma^{2}_{\delta_{\omega}}\mathcal{I}_{3}) .

RMSE		\adl@mkpreamc\@addtopreamble   \@arstrut\@preamble				\adl@mkpreamc\@addtopreamble      \@arstrut\@preamble
	$x$	$y$	$z$		$x$	$y$	$z$
ML	5.0	1.0	-4.0		5.1	5.3	6.4
MAP ${\sigma }_{{\delta }_{\omega }}=0.05$	5.0	1.0	-4.0		5.1	5.3	6.4
MAP ${\sigma }_{{\delta }_{\omega }}=1\cdot {10}^{-3}$	3.9	0.8	-2.8		4.1	4.0	4.7

Equations

$$x^{\text{n}}=R^{\text{nb}}{x}^{\text{b}},x^{\text{b}}=(R^{\text{nb}}{)}^{\mathsf{T}}{x}^{\text{n}}=R^{\text{bn}}{x}^{\text{n}}.$$

$${\omega }_{\text{ib}}^{\text{b}}=R^{\text{bn}}\left({\omega }_{\text{ie}}^{\text{n}}+{\omega }_{\text{en}}^{\text{n}}\right)+{\omega }_{\text{nb}}^{\text{b}},$$

$$f^{\text{b}}=R^{\text{bn}}(a_{\text{ii}}^{\text{n}}-g^{\text{n}}),$$

$$a_{\text{ii}}^{\text{n}}=R^{\text{ne}}{R}^{\text{ei}}{a}_{\text{ii}}^{\text{i}}.$$

$${\frac{\text{d}}{\text{d}t}{p}^{\text{n}}|}_{\text{n}}=v_{\text{n}}^{\text{n}},{\frac{\text{d}}{\text{d}t}{v}^{\text{n}}|}_{\text{n}}=a_{\text{nn}}^{\text{n}}.$$

$${\frac{\text{d}}{\text{d}t}{x}^{\text{u}}|}_{\text{u}}={\frac{\text{d}}{\text{d}t}{R}^{\text{uv}}{x}^{\text{v}}|}_{\text{u}}=R^{\text{uv}}{\frac{\text{d}}{\text{d}t}{x}^{\text{v}}|}_{\text{v}}+{\omega }_{\text{uv}}^{\text{u}}\times x^{\text{u}},$$

$$p^{\text{i}}=R^{\text{ie}}{p}^{\text{e}},$$

$$v_{\text{i}}^{\text{i}}$$

$$a_{\text{ii}}^{\text{i}}$$

$$=a_{\text{ee}}^{\text{i}}+2{\omega }_{\text{ie}}^{\text{i}}\times v_{\text{e}}^{\text{i}}+{\omega }_{\text{ie}}^{\text{i}}\times {\omega }_{\text{ie}}^{\text{i}}\times p^{\text{i}},$$

$$p^{\text{e}}=R^{\text{en}}{p}^{\text{n}}+n_{\text{ne}}^{\text{e}},$$

$$v_{\text{e}}^{\text{e}}$$

$$a_{\text{ee}}^{\text{e}}$$

$$a_{\text{ii}}^{\text{n}}=a_{\text{nn}}^{\text{n}}+2{\omega }_{\text{ie}}^{\text{n}}\times v_{\text{n}}^{\text{n}}+{\omega }_{\text{ie}}^{\text{n}}\times {\omega }_{\text{ie}}^{\text{n}}\times p^{\text{n}},$$

$$\Vert {\omega }_{\text{ie}}^{\text{n}}\times {\omega }_{\text{ie}}^{\text{n}}\times p^{\text{n}}{\Vert }_2\le \Vert {\omega }_{\text{ie}}^{\text{n}}{\Vert }_2\Vert {\omega }_{\text{ie}}^{\text{n}}{\Vert }_2\Vert p^{\text{n}}{\Vert }_2.$$

$$y_t^1={\left(\begin{array}{cc}0& 0\end{array}\right)}^{\mathsf{T}},y_t^2={\left(\begin{array}{cc}2& 0\end{array}\right)}^{\mathsf{T}}.$$

$$y_t^1$$

$$y_t^2$$

$$p_t\sim \mathcal{N}\left(\left(\begin{array}{c}0.4\\ 0\end{array}\right),0.2{\mathcal{I}}_2\right).$$

$$p(x_{1:N},\theta \mid y_{1:N}),$$

$$p(x_t\mid y_{1:t}).$$

$$p(x_{1:N},\theta \mid y_{1:N})$$

$$p(x_t\mid y_{1:t})$$

$$p(x_t\mid y_{1:t-1})=\int p(x_t\mid x_{t-1})p(x_{t-1}\mid y_{1:t-1})\mathsf{d}{x}_{t-1}.$$

$$x_{t+1}=f_t(x_t,w_t).$$

$$p(x_{t+1}\mid x_t)\sim \mathcal{N}(x_{t+1};f_t(x_t),Q),$$

$$y_t=h_t(x_t,e_t),$$

$$R{R}^{\mathsf{T}}=R^{\mathsf{T}}R={\mathcal{I}}_3,\det R=1.$$

$$x^{\text{u}}$$

$$x^{\text{v}}$$

$$x_{\parallel }^{\text{v}}$$

$$x_{\star }^{\text{v}}$$

$${\left(x_{\star }^{\text{v}}\right)}_{\parallel }$$

Full text

1. 1 Introduction

2. 1.1 Background and motivation

3. 1.2 Using inertial sensors for position and orientation estimation

4. 1.3 Tutorial content and its outline

5. 2 Inertial Sensors

6. 2.1 Coordinate frames

7. 2.2 Angular velocity

8. 2.3 Specific force

9. 2.4 Sensor errors

10. 3 Probabilistic Models

11. 3.1 Introduction

12. 3.2 Parametrizing orientation

13. 3.2.1 Rotation matrices

14. 3.2.2 Rotation vector

15. 3.2.3 Euler angles

16. 3.2.4 Unit quaternions

17. 3.3 Probabilistic orientation modeling

18. 3.3.1 Linearization

19. 3.3.2 Alternative methods

20. 3.4 Measurement models

21. 3.4.1 Gyroscope measurement models

22. 3.4.2 Accelerometer measurement models

23. 3.4.3 Modeling additional information

24. 3.5 Choosing the state and modeling its dynamics

25. 3.6 Models for the prior

26. 3.7 Resulting probabilistic models

27. 3.7.1 Pose estimation

28. 3.7.2 Orientation estimation

29. 4 Estimating Position and Orientation

30. 4.1 Smoothing in an optimization framework

31. 4.1.1 Gauss-Newton optimization

32. 4.1.2 Smoothing estimates of the orientation using optimization

33. 4.1.3 Computing the uncertainty

34. 4.2 Filtering estimation in an optimization framework

35. 4.2.1 Filtering estimates of the orientation using optimization

36. 4.3 Extended Kalman filtering

37. 4.3.1 Estimating orientation using quaternions as states

38. 4.3.2 Estimating orientation using orientation deviations as states

39. 4.4 Complementary filtering

40. 4.5 Evaluation based on experimental and simulated data

41. 4.5.1 General characteristics

42. 4.5.2 Representing uncertainty

43. 4.5.3 Comparing the different algorithms

44. 4.6 Extending to pose estimation

45. 5 Calibration

46. 5.1 Maximum a posteriori calibration

47. 5.2 Maximum likelihood calibration

48. 5.3 Orientation estimation with an unknown gyroscope bias

49. 5.4 Identifiability

50. 6 Sensor Fusion Involving Inertial Sensors

51. 6.1 Non-Gaussian noise distributions: time of arrival measurements

52. 6.2 Using more complex sensor information: inertial and vision

53. 6.3 Including non-sensory information: inertial motion capture

54. 6.4 Beyond pose estimation: activity recognition

55. 7 Concluding Remarks

56. 8 Notation and Acronyms

57. A Orientation Parametrizations

58. A.1 Quaternion algebra

59. A.2 Conversions between different parametrizations

60. B Pose Estimation

61. B.1 Smoothing in an optimization framework

62. B.2 Filtering in an optimization framework

63. B.3 Extended Kalman filter with quaternion states

64. B.4 Extended Kalman filter with orientation deviation states

65. C Gyroscope Bias Estimation

66. C.1 Smoothing in an optimization framework

67. C.2 Filtering in an optimization framework

68. C.3 Extended Kalman filter with quaternion states

69. C.4 Extended Kalman filter with orientation deviation states

Chapter 1 Introduction

In this tutorial, we discuss the topic of position and orientation estimation using inertial sensors. We consider two separate problem formulations. The first is estimation of orientation only, while the other is the combined estimation of both position and orientation. The latter is sometimes called pose estimation. We start by providing a brief background and motivation in §1.1 and explain what inertial sensors are and give a few concrete examples of relevant application areas of pose estimation using inertial sensors. In §1.2, we subsequently discuss how inertial sensors can be used to provide position and orientation information. Finally, in §1.3 we provide an overview of the contents of this tutorial as well as an outline of subsequent chapters.

1.1 Background and motivation

The term inertial sensor is used to denote the combination of a three-axis accelerometer and a three-axis gyroscope. Devices containing these sensors are commonly referred to as inertial measurement units (IMUs). Inertial sensors are nowadays also present in most modern smartphone, and in devices such as Wii controllers and virtual reality (VR) headsets, as shown in Figure LABEL:fig:intro-inertialSensors.

A gyroscope measures the sensor’s angular velocity, i.e. the rate of change of the sensor’s orientation. An accelerometer measures the external specific force acting on the sensor. The specific force consists of both the sensor’s acceleration and the earth’s gravity. Nowadays, many gyroscopes and accelerometers are based on microelectromechanical system (MEMS) technology. MEMS components are small, light, inexpensive, have low power consumption and short start-up times. Their accuracy has significantly increased over the years.

There is a large and ever-growing number of application areas for inertial sensors, see e.g. [7, 59, 109, 156]. Generally speaking, inertial sensors can be used to provide information about the pose of any object that they are rigidly attached to. It is also possible to combine multiple inertial sensors to obtain information about the pose of separate connected objects. Hence, inertial sensors can be used to track human motion as illustrated in Figure 1.2. This is often referred to as motion capture. The application areas are as diverse as robotics, biomechanical analysis and motion capture for the movie and gaming industries. In fact, the use of inertial sensors for pose estimation is now common practice in for instance robotics and human motion tracking, see e.g. [86, 54, 112]. A recent survey [1] shows that 28% of the contributions to the IEEE International Conference on Indoor Positioning and Indoor Navigation (IPIN) make use of inertial sensors. Inertial sensors are also frequently used for pose estimation of cars, boats, trains and aerial vehicles, see e.g. [139, 23]. Examples of this are shown in Figure 1.3.

There exists a large amount of literature on the use of inertial sensors for position and orientation estimation. The reason for this is not only the large number of application areas. Important reasons are also that the estimation problems are nonlinear and that different parametrizations of the orientation need to be considered [47, 77], each with its own specific properties. Interestingly, approximative and relatively simple position and orientation estimation algorithms work quite well in practice. However, careful modeling and a careful choice of algorithms do improve the accuracy of the estimates.

In this tutorial we focus on the signal processing aspects of position and orientation estimation using inertial sensors, discussing different modeling choices and a number of important algorithms. These algorithms will provide the reader with a starting point to implement their own position and orientation estimation algorithms.

1.2 Using inertial sensors for position and orientation estimation

As illustrated in §1.1, inertial sensors are frequently used for navigation purposes where the position and the orientation of a device are of interest. Integration of the gyroscope measurements provides information about the orientation of the sensor. After subtraction of the earth’s gravity, double integration of the accelerometer measurements provides information about the sensor’s position. To be able to subtract the earth’s gravity, the orientation of the sensor needs to be known. Hence, estimation of the sensor’s position and orientation are inherently linked when it comes to inertial sensors. The process of integrating the measurements from inertial sensors to obtain position and orientation information, often called dead-reckoning, is summarized in Figure 1.4.

If the initial pose would be known, and if perfect models for the inertial sensor measurements would exist, the process illustrated in Figure 1.4 would lead to perfect pose estimates. In practice, however, the inertial measurements are noisy and biased as will be discussed in more detail in §2.4. Because of this, the integration steps from angular velocity to rotation and from acceleration to position introduce integration drift. This is illustrated in Example 1.1.

Example 1.1** (Integration drift)**

Let us first focus on the general case of measuring a quantity that is constant and equal to zero. The integrated and double integrated signals are therefore also equal to zero. However, let us now assume that we measure this quantity using a non-perfect sensor. In case our measurements are corrupted by a constant bias, integration of these measurements will lead to a signal which grows linearly with time. Double integration leads to a signal that instead grows quadratically with time. If the sensor instead measures a zero-mean white noise signal, the expected value of the integrated measurements would be zero, but the variance would grow with time. This is illustrated in Figure 1.5 for the integration of a signal $y_{t}=e_{t}$ with $e_{t}\sim\mathcal{N}(0,1)$. Hence, integration drift is both due to integration of a constant bias and due to integration of noise.

To illustrate integration drift using experimental data, a stationary data set is collected with a Sony Xperia Z5 Compact smartphone using the app described in [57]. The smartphone contains accelerometers and gyroscopes produced by Invensense [64]. We integrate the inertial measurements to obtain position and orientation estimates. Since the smartphone is kept stationary during the data collection, we expect the position and orientation to remain the same. However, the orientation estimates drift a few degrees over $10$ seconds as shown in Figure 1.6(a). Note that the integration drift is not the same for all axes. This is mainly due to a different sensor bias in the different axes. This will be studied further in Example 2.3, where the same data set is used to study the sensor characteristics. As shown in Figure 1.6(b), the position drifts several meters over $10\mathrm{s}$. The reason for this is two-fold. First, the accelerometer measurements need to be integrated twice. Second, the orientation estimates need to be used to subtract the gravity and any errors in this will result in leakage of gravity into the other components.

**

From the example above, it can be concluded that errors in the measurements have a large impact on the quality of the estimated position and orientation using inertial sensors only. This is particularly the case for position, which relies both on double integration of the acceleration and on accurate orientation estimates to subtract the earth’s gravity. Because of this, inertial sensors need to be supplemented with other sensors and other models to obtain accurate position and orientation estimates.

Inertial sensors provide pose estimates at high sampling rates which are accurate on a short time scale but drift over longer time scales. They are therefore very suitable for being combined with sensors with a lower sampling rate but with information that does not drift over time. For pose estimation, inertial sensors are often combined with measurements from for instance a global navigation satellite system (GNSS) [71, 147, 59], an ultrawideband (UWB) system [74, 132, 110, 27, 30] or cameras [26, 61, 79, 94]. For orientation estimation, they are often used in combination with magnetometers, which measure the direction of the magnetic field [124, 120].

This tutorial aims at giving an introduction on how to use inertial sensors for position and orientation estimation, but also on how to combine them with additional information. These additional sensors are not the focus of this paper but simple models will be used for magnetometers and sensors providing position information to illustrate the combined use of these sensors.

1.3 Tutorial content and its outline

To obtain accurate position and orientation estimates using inertial sensors in combination with additional measurements and models, a number of important things need to be considered. First, the quantities measured by the inertial sensors need to be accurately described and the sources of error need to be characterized. This is the topic of Chapter 2. Note that throughout the tutorial, we will focus on MEMS inertial sensors and consider both data from standalone IMUs and from smartphones. This implies that we do not focus on for instance mechanical or optical gyroscopes and on mechanical or solid-state accelerometers [147]. These sensors may have characteristics that are quite different from the MEMS inertial sensors considered here.

Based on the analysis of the sensors in Chapter 2 and on additional analysis of the application at hand, models can be constructed. This is the topic of Chapter 3, where we will also discuss different parametrizations of orientation. This will highlight the challenges in parametrizing and estimating orientations and show that the orientation estimation problem is inherently nonlinear. Furthermore, we will present two models that can be used for position and orientation estimation. The first is a model for pose estimation using inertial measurements in combination with position measurements. The second is a model for orientation estimation, using inertial and magnetometer measurements.

In Chapter 4, different algorithms for position and orientation estimation will be introduced. The general structure of the algorithms will be discussed, after which explicit algorithms for orientation estimation using inertial and magnetometer measurements are given. We will also discuss how the algorithms can be extended to pose estimation when position measurements are available. Some general characteristics of the two estimation problems will be given and the quality of the estimates from the different algorithms will be analyzed. Which algorithm is most suitable for which application depends strongly on the computational power that is available, the accuracy that is required and the characteristics of the problem at hand.

In Chapter 4, we assume that the sensors are properly calibrated. However, calibration of the sensors is important to for instance estimate the inertial sensor biases. Furthermore, calibration is specifically of concern when combining inertial data with other sensors. In these cases, it is important that the inertial sensor axes and the axes of the additional sensors are aligned. Sensor calibration is the topic of Chapter 5. As an illustrative example, we will consider the estimation of an unknown gyroscope bias.

Our focus in this tutorial is to present models and algorithms that can be used for position and orientation estimation using inertial measurements. Because of this, we assume fairly simple models for the additional information — the magnetometer and position measurements. In Chapter 6, we will discuss how the algorithms from Chapter 4 can be used in more complex settings. For example, we will consider the cases of more complex position information in the form of images from a camera, the presence of non-Gaussian measurement noise and the availability of additional information that can be exploited. The information provided by the inertial sensors remains one of the main building blocks of algorithms that can be used for these cases. Adaptations of the algorithms presented in Chapter 4 can therefore be used to also solve these more complex scenarios. We will end this tutorial with some concluding remarks in Chapter 7.

Chapter 2 Inertial Sensors

To combine inertial measurements with additional sensors and models for position and orientation estimation, it is important to accurately describe the quantities measured by the inertial sensors as well as to characterize the typical sensor errors. This will be the topic of this chapter. It will serve as a basis for the probabilistic models discussed in Chapter 3.

As discussed in Chapter 1, accelerometers and gyroscopes measure the specific force and the angular velocity, respectively. In §2.2 and §2.3, we will discuss these quantities in more detail. To enable a discussion about this, in §2.1 a number of coordinate frames and the transformations between them will be discussed. We assume that we have 3D accelerometers and 3D gyroscopes, i.e. that the sensors have three sensitive axes along which these physical quantities are measured. They are measured in terms of an output voltage which is converted to a physical measurement based on calibration values obtained in the factory. Even though the sensors are typically calibrated in the factory, (possibly time-varying) errors can still remain. In §2.4, the most commonly occurring sensor errors are discussed.

2.1 Coordinate frames

In order to discuss the quantities measured by the accelerometer and gyroscope in more detail, a number of coordinate frames need to be introduced:

The body frame $\boldsymbol{b}$

is the coordinate frame of the moving IMU. Its origin is located in the center of the accelerometer triad and it is aligned to the casing. All the inertial measurements are resolved in this frame.

The navigation frame $\boldsymbol{n}$

is a local geographic frame in which we want to navigate. In other words, we are interested in the position and orientation of the $b$-frame with respect to this frame. For most applications it is defined stationary with respect to the earth. However, in cases when the sensor is expected to move over large distances, it is customary to move and rotate the $n$-frame along the surface of the earth. The first definition is used throughout this tutorial, unless mentioned explicitly.

The inertial frame $\boldsymbol{i}$

is a stationary frame. The IMU measures linear acceleration and angular velocity with respect to this frame. Its origin is located at the center of the earth and its axes are aligned with respect to the stars.

The earth frame $\boldsymbol{e}$

coincides with the $i$-frame, but rotates with the earth. That is, it has its origin at the center of the earth and axes which are fixed with respect to the earth.

The $n$, $i$ and $e$ coordinate frames are illustrated in Figure 2.1. We use a superscript to indicate in which coordinate frame a vector is expressed. Vectors can be rotated from one coordinate frame to another using a rotation matrix. We use a double superscript to indicate from which coordinate frame to which coordinate frame the rotation is defined. An illustration is given in Example 2.1.

Example 2.1** (Rotation of vectors to different coordinate frames)**

Consider a vector $x$ expressed in the body frame $b$. We denote this vector by $x^{\text{b}}$. The rotation matrix $R^{\text{nb}}$ rotates a vector from the body frame $b$ to the navigation frame $n$. Conversely, the rotation from navigation frame $n$ to body frame $b$ is denoted $R^{\text{bn}}=(R^{\text{nb}})^{\mathsf{T}}$. Hence, the vector $x$ expressed in the body frame ($x^{\text{b}}$) and expressed in the navigation frame ($x^{\text{n}}$) are related according to

[TABLE]

**

2.2 Angular velocity

The gyroscope measures the angular velocity of the body frame with respect to the inertial frame, expressed in the body frame [147], denoted by $\omega_{\text{ib}}^{\text{b}}$. This angular velocity can be expressed as

[TABLE]

where $R^{\text{bn}}$ is the rotation matrix from the navigation frame to the body frame. The earth rate, i.e. the angular velocity of the earth frame with respect to the inertial frame is denoted by $\omega_{\text{ie}}$. The earth rotates around its own $z$-axis in $23.9345$ hours with respect to the stars [101]. Hence, the earth rate is approximately $7.29\cdot 10^{-5}~{}\radianpersecond$.

In case the navigation frame is not defined stationary with respect to the earth, the angular velocity $\omega_{\text{en}}$, i.e. the transport rate is non-zero. The angular velocity required for navigation purposes — in which we are interested when determining the orientation of the body frame with respect to the navigation frame — is denoted by $\omega_{\text{nb}}$.

2.3 Specific force

The accelerometer measures the specific force $f$ in the body frame $b$ [147]. This can be expressed as

[TABLE]

where $g$ denotes the gravity vector and $a_{\text{ii}}^{\text{n}}$ denotes the linear acceleration of the sensor expressed in the navigation frame, which is

[TABLE]

The subscripts on the linear acceleration $a$ are used to indicate in which frame the differentiation is performed. For navigation purposes, we are interested in the position of the sensor in the navigation frame $p^{\text{n}}$ and its derivatives as performed in the navigation frame

[TABLE]

A relation between $a_{\text{ii}}$ and $a_{\text{nn}}$ can be derived by using the relation between two rotating coordinate frames. Given a vector $x$ in a coordinate frame $u$,

[TABLE]

where $\omega_{\text{uv}}^{\text{u}}$ is the angular velocity of the $v$-frame with respect to the $u$-frame, expressed in the $u$-frame. For a derivation of this relation in the context of inertial navigation, see [59, 147]. For a general introduction, see any textbook on dynamics, e.g. [92, 96].

Using the fact that

[TABLE]

the velocity $v_{\text{i}}$ and acceleration $a_{\text{ii}}$ can be expressed as

[TABLE]

where we have made use of (2.5), (2.6), and the fact that the angular velocity of the earth is constant, i.e. $\tfrac{\textrm{d}}{\textrm{d}t}\omega_{\text{ie}}^{\text{i}}=0$. Using the relation between the earth frame and the navigation frame

[TABLE]

where $n_{\text{ne}}$ is the distance from the origin of the earth coordinate frame to the origin of the navigation coordinate frame, expressions similar to (2.8) can be derived. Note that in general it can not be assumed that $\tfrac{\textrm{d}}{\textrm{d}t}\omega_{\text{en}}=0$. Inserting the obtained expressions into (2.8), it is possible to derive the relation between $a_{\text{ii}}$ and $a_{\text{nn}}$. Instead of deriving these relations, we will assume that the navigation frame is fixed to the earth frame, and hence $R^{\text{en}}$ and $n_{\text{ne}}^{\text{e}}$ are constant and

[TABLE]

This is a reasonable assumption as long as the sensor does not travel over significant distances as compared to the size of the earth and it will be one of the model assumptions that we will use in this tutorial. More on the modeling choices will be discussed in Chapter 3.

Inserting (2.10) into (2.8) and rotating the result, it is possible to express $a_{\text{ii}}^{\text{n}}$ in terms of $a_{\text{nn}}^{\text{n}}$ as

[TABLE]

where $a_{\text{nn}}$ is the acceleration required for navigation purposes. The term $\omega_{\text{ie}}^{\text{n}}\times\omega_{\text{ie}}^{\text{n}}\times p^{\text{n}}$ is known as the centrifugal acceleration and $2\omega_{\text{ie}}^{\text{n}}\times v^{\text{n}}_{\text{n}}$ is known as the Coriolis acceleration. The centrifugal acceleration is typically absorbed in the (local) gravity vector. In Example 2.2, we illustrate the magnitude of both the centrifugal and the Coriolis acceleration.

Example 2.2** (Magnitude of centrifugal and Coriolis acceleration)**

The centrifugal acceleration depends on the location on the earth. It is possible to get a feeling for its magnitude by considering the property of the cross product stating that

[TABLE]

Since the magnitude of $\omega_{\text{ie}}$ is approximately $7.29\cdot 10^{-5}~{}\mathrm{rad}\mathrm{/}\mathrm{s}$ and the average radius of the earth is $6371~{}\mathrm{k}\mathrm{m}$ [101], the magnitude of the centrifugal acceleration is less than or equal to $3.39\cdot 10^{-2}~{}\metrepersquaresecond$.

The Coriolis acceleration depends on the speed of the sensor. Let us consider a person walking at a speed of $5~{}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{h}$. In that case the magnitude of the Coriolis acceleration is approximately $2.03\cdot 10^{-4}~{}\metrepersquaresecond$. For a car traveling at $120~{}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{h}$, the magnitude of the Coriolis acceleration is instead $4.86\cdot 10^{-3}~{}\metrepersquaresecond$.

**

2.4 Sensor errors

As discussed in §2.2 and §2.3, the gyroscope measures the angular velocity $\omega^{\text{b}}_{\text{ib}}$ and the accelerometer measures the specific force $f^{\text{b}}$. However, as already briefly mentioned in §1.2, there are several reasons for why this is not exactly the case. Two of these reasons are a slowly time-varying sensor bias and the presence of measurement noise. The sensor errors in the inertial measurements are illustrated in Example 2.3 using experimental data.

Example 2.3** (Inertial sensor measurements and their errors)**

In Figures 2.2–2.4, gyroscope and accelerometer measurements are displayed for around $10$ seconds of stationary data collected with a Sony Xperia Z5 Compact smartphone. Since the smartphone is stationary, the gyroscope is expected to only measure the earth’s angular velocity. However, as can be seen in Figure 2.2(a), the gyroscope measurements are corrupted by noise. As shown in Figure 2.3, this noise can be seen to be quite Gaussian. Furthermore, the measurements can be seen to be biased.

During the stationary period, we would expect the accelerometer to measure the gravity, the centrifugal acceleration and the Coriolis acceleration. Note that again the measurements are corrupted by noise, which can be seen to be quite Gaussian in Figure 2.4. The $x$- and $y$-components of the accelerometer measurements are not zero-mean. This can be due to the fact that the table on which the smartphone lies is not completely flat, implying that part of the gravity vector is measured in these components. It can also reflect a sensor bias. The $z$-component is actually larger than expected which indicates the presence of an accelerometer bias at least in this axis.

Note that from the above discussion it can be concluded that it is more straightforward to determine the gyroscope bias than it is to determine the accelerometer bias. To be able to estimate the gyroscope bias, it is sufficient to leave the sensor stationary. For the accelerometer, however, it is in that case difficult to distinguish between a bias and a table that is not completely flat. For more information about accelerometer or general inertial sensor calibration, see [145, 106, 108].

The gyroscope in the smartphone is automatically recalibrated during stationary time periods. The measurements shown in Figure 2.2(a) have not been corrected for this (so-called uncalibrated or raw data). The reason why the gyroscope is calibrated during stationary periods, is because its bias is slowly time-varying. As an example, let us consider a data set of $55$ minutes. The gyroscope bias during the first minute of the data set was $\begin{pmatrix}35.67&56.22&0.30\end{pmatrix}^{\mathsf{T}}\cdot 10^{-4}~{}\radianpersecond$, while the gyroscope bias during the last minute of the data set was $\begin{pmatrix}37.01&53.17&-1.57\end{pmatrix}^{\mathsf{T}}\cdot 10^{-4}~{}\radianpersecond$.

**

The performance of IMUs is often specified in terms of their Allan variance [63, 33, 2]. The Allan variance gives information about the sensor errors for stationary conditions, i.e. in a stable climate without exciting the system. It studies the effect of averaging measurements for different cluster times $T_{c}$. Typically, the Allan standard deviation $\sigma_{\text{A}}(T_{c})$ is plotted against the cluster time $T_{c}$ as illustrated in Figure 2.5. This figure shows the characteristic behavior of the Allan variance for inertial sensors. To study it more in detail, we will discuss two components of the Allan variance that are typically of concern for inertial sensors: the white noise and the bias instability.

Assume, as in Example 1.1, that we have a white noise signal with standard deviation $\sigma$. A longer averaging time would for this signal lead to values closer to zero. The contribution to the Allan standard deviation from the white noise component is given by $\sigma_{\text{A}}(T_{c})=\tfrac{\sigma}{\sqrt{n}}$ where $n$ is the number of samples averaged over. This corresponds to a line with slope $-1/2$ in a $\log$–$\log$ plot. For instance in the Allan deviation for the gyroscope in Figure 2.5, the lines can be seen to have a slope of $-1/2$ until around $10-20~{}\mathrm{s}$, which indicates that the white noise is the dominating source of error for these short integration times.

A constant bias does not have any effect on the Allan variance diagram. However, in case the bias changes, longer averaging times will no longer be beneficial. Hence, the Allan variance diagrams in Figure 2.5 show a deviation from the slope $-1/2$ for longer averaging times.

The Allan variance is a useful tool to study and compare the noise characteristics of inertial sensors. However, it only considers stationary conditions. In dynamic conditions, a large number of other error sources potentially come into play, see e.g. [147, 154]. These are for instance related to the fact that the sensors sample at discrete times. Hence, to capture high-frequency signals, high sampling frequencies are desired [128, 129]. Furthermore, large dynamics can lead to erroneous or saturated measurements. Other errors that are not included are for instance changes in the sensitivity of the axes due to changes in temperature. We should therefore never just rely on the Allan variance when deciding which sensor to use in a particular application.

Chapter 3 Probabilistic Models

Pose estimation is about estimating the position and orientation of the body frame $b$ in the navigation frame $n$. This problem is illustrated in Figure 3.1, where the position and orientation of the body changes from time $t_{1}$ to time $t_{2}$. In this chapter, we will introduce the concept of probabilistic models and discuss different modeling choices when using inertial sensors for pose estimation.

The subject of probabilistic modeling is introduced in §3.1. Most complexity in pose estimation lies in the nonlinear nature of the orientation and the fact that orientation can be parametrized in different ways. How to parametrize the orientation is therefore a crucial modeling choice in any pose estimation algorithm. Because of this, we will discuss different parametrizations for the orientation in §3.2 and in §3.3 we will discuss how these different parametrizations can be used in probabilistic modeling.

Our probabilistic models consist of three main components. First, in §3.4, we introduce models describing the knowledge about the pose that can be inferred from the measurements. Second, in §3.5, we model how the sensor pose changes over time. Finally, in §3.6, models of the initial pose are introduced.

The chapter will conclude with a discussion on the resulting probabilistic models in §3.7. The models that will be used in the position and orientation estimation algorithms in Chapter 4 will be introduced in this section.

3.1 Introduction

Probabilistic models constitute the foundation of the estimation algorithms in Chapter 4. In this section we will introduce the concept of probabilistic modeling and the notation that is used in building our models. Models are used to describe the information about the dynamics and the available measurements. These models are subsequently used in combination with the measurements to infer some knowledge. The knowledge that we are interested in is the pose of the sensor and we use information about the sensor dynamics and the available measurements (amongst others, inertial measurements). A simplified case where probabilistic modeling is used to estimate the position of a sensor is given in Example 3.1.

Example 3.1** (Probabilistic modeling)**

Let us estimate the 2D position $p_{t}$ of a sensor at time $t$ from two position measurements

[TABLE]

A straightforward suggestion for an estimate of the position would be $\hat{p}_{t}=\begin{pmatrix}1&0\end{pmatrix}^{\mathsf{T}}$. Let us now assume that we know the accuracy of the sensors and represent this in terms of the following probabilistic models

[TABLE]

where $\mathcal{I}_{2}$ denotes a $2\times 2$ identity matrix. A reasonable position estimate would instead be

[TABLE]

We will not go into details about how this estimate is derived. Instead, we would like to point out two differences between this position estimate and our initial suggestion $\hat{p}_{t}$. First, based on the knowledge of the accuracy of the sensors, it is sensible to trust the measurement from the first sensor more than the measurement from the second sensor. Our improved position estimate is therefore closer to the measurement of the first sensor than our initial suggestion $\hat{p}_{t}$. Furthermore, based on the accuracy of the sensors, it is possible to derive the accuracy of our estimate.

Now consider the case where we are also interested in estimating the position $p_{t+1}$. Knowledge that the sensor is worn by a human or placed in a car, would give us information about how far the sensor can travel from time $t$ to time $t+1$. If the sensor would be placed in for instance a train, the motion would even be constrained to be along the tracks. Incorporating this information about the dynamics of the sensor will improve the estimate of $p_{t+1}$.

**

We split the knowledge that we want to infer into the unknown time-varying states $x_{t}$ for $t=1,\ldots,N$, or equivalently $x_{1:N}$, and the unknown constant parameters $\theta$. We denote the measurements by $y_{k}$ for $k=1,\ldots,K$. The times $k$ at which these measurements are obtained do not necessarily correspond with the times $t$ at which the states are defined. It is also not necessary for all sensors to sample at the same frequency. As discussed in §2.4, the inertial sensors are typically sampled at fairly high rates to capture high-frequency dynamics. In stand-alone, wired IMUs, all sensors typically have the same, constant sampling frequency. Specifically in the case of wireless sensors and smartphones, however, the sampling frequencies can vary both over sensors and over time. In the remainder, we assume that the times $t$ at which the states are defined coincide with the times $k$ at which the gyroscopes sample. Hence, we denote the gyroscope measurements $y_{\omega,t}$ with $t=1,\ldots N$. For notational convenience, we will also use the subscript $t$ for the measurements from other sensors. Note that these are not required to actually sample at each time $t$ for $t=1,\ldots N$. For instance, magnetometers in smartphones often sample either at equal or at half the sampling frequencies of the inertial sensors, while position aiding sensors like for instance GNSS or UWB typically sample at much lower sampling frequencies.

Our aim is now to infer information about the states $x_{1:N}$ and the parameters $\theta$ using the measurements $y_{1:N}$ and the probabilistic models. This can be expressed in terms of a conditional probability distribution

[TABLE]

where $p(a\mid b)$ denotes the conditional probability of $a$ given $b$. In the pose estimation problem, we are interested in obtaining point estimates which we denote $\hat{x}_{1:N}$ and $\hat{\theta}$. It is typically also highly relevant to know how certain we are about these estimates. This is often expressed in terms of a covariance. When the distribution (3.1) is Gaussian, the distribution is completely described in terms of its mean and covariance.

In (3.1) we assume that all measurements $y_{1:N}$ are used to obtain the posterior distribution of $x_{1:N}$ and $\theta$. This is referred to as smoothing. Although it makes sense to use all available information to obtain the best estimates, a downside of smoothing is that we need to wait until all measurements are collected before the pose can be computed. Because of this, in many applications, we are also interested in filtering. In filtering we estimate $x_{t}$ using all measurements up to and including time $t$. One way of dealing with constant parameters in filtering is to treat them as slowly time-varying. In this case, they can be considered to be included in the time-varying states $x_{t}$. The filtering problem can be expressed in terms of the conditional probability distribution

[TABLE]

We have now introduced smoothing, where the states $x_{1:N}$ are estimated simultaneously, and filtering, where at each time instance the state $x_{t}$ is estimated. There is a large range of intermediate methods, where a batch of states $x_{t-L_{1}:t+L_{2}}$, with $L_{1}$ and $L_{2}$ being positive integers, is estimated using the measurements $y_{1:t}$. This is related to fixed-lag smoothing and moving horizon estimation [66, 113].

The topic of how to estimate the conditional probability distributions for position and orientation estimation will be introduced in Chapter 4. We will now instead take a closer look at these distributions and their different components. A fundamental assumption here is that we assume that our models possess the Markov property, implying that all information up to the current time $t$ is contained in the state $x_{t}$. This is illustrated in Figure 3.2 in terms of a probabilistic graphical model [12]. The state $x_{t+1}$ can be seen to depend on $x_{t}$ and to result in the measurements $y_{t+1}$. It is conditionally independent of $x_{1:t-1}$ given the state $x_{t}$. Using Bayes’ rule and the Markov property, the conditional distributions (3.1) and (3.2) can be decomposed as

[TABLE]

The predictive distribution $p(x_{t}\mid y_{1:t-1})$ can be computed by marginalizing out the previous state $x_{t-1}$ as

[TABLE]

In (3.3), $p(\theta)$ and $p(x_{1}\mid\theta)$ encode our prior information of $\theta$ and the knowledge of the state $x_{1}$ given $\theta$, respectively. The dynamics are modeled in terms of $p(x_{t+1}\mid x_{t},\theta)$ and $p(x_{t+1}\mid x_{t})$. The distributions $p(y_{t}\mid x_{t},\theta)$ and $p(y_{t}\mid x_{t})$ model the information given by the measurements about the state and the parameters.

The dynamics of the state can be modeled in terms of a nonlinear function $f_{t}(\cdot)$ as

[TABLE]

The uncertainty of the dynamic model is modeled in terms of $w_{t}$, which is often referred to as the process noise. The model (3.5) provides information about the distribution $p(x_{t+1}\mid x_{t})$. More explicitly, if $w_{t}$ is Gaussian additive noise with $w_{t}\sim\mathcal{N}(0,Q)$, then

[TABLE]

where we use the notation $\mathcal{N}(x_{t+1}\,;\,f_{t}(x_{t}),Q)$ to explain that the random variable $x_{t+1}$ is normal distributed with mean $f_{t}(x_{t})$ and covariance $Q$.

The information given by the measurements about the state $x_{t}$ can be modeled as

[TABLE]

where $h_{t}(\cdot)$ is a possibly nonlinear function and $e_{t}$ is the measurement noise. The measurement model (3.7) provides information about the distribution $p(y_{t}\mid x_{t})$. The combination of (3.5), (3.7) and a model of the prior $p(x_{1})$ is referred to as a state space model [69] which is widely used in a large number of fields.

3.2 Parametrizing orientation

Rotating a vector in $\mathbb{R}^{3}$ changes the direction of the vector while retaining its length. The group of rotations in $\mathbb{R}^{3}$ is the special orthogonal group $SO(3)$. In this section we introduce four different ways of parametrizing orientations. Note that these describe the same quantity and can hence be used interchangeably. The different parametrizations can be converted to one another, see also Appendix A. There are differences in for instance the number of parameters used in the representation, the singularities and the uniqueness.

3.2.1 Rotation matrices

We encountered rotation matrices already in Chapter 2. Rotation matrices $R\in\mathbb{R}^{3\times 3}$ have the following properties

[TABLE]

The properties (3.8) provide an interpretation of the name special orthogonal group $SO(3)$. All orthogonal matrices of dimension $3\times 3$ have the property $RR^{\mathsf{T}}=R^{\mathsf{T}}R=\mathcal{I}_{3}$ and are part of the orthogonal group $O(3)$. The notion special in $SO(3)$ specifies that only matrices with $\det R=1$ are considered rotations.

Consider two coordinate frames denoted $u$ and $v$. As was illustrated in Example 2.1, a vector $x$ expressed in the $v$-frame can be rotated to the $u$-frame as

[TABLE]

A rotation matrix is a unique description of the orientation. It has $9$ components which depend on each other as defined in (3.8).

3.2.2 Rotation vector

As described by Leonhard Euler in [34], a rotation around a point is always equivalent to a single rotation around some axis through this point, see [107] for a number of proofs. This is generally referred to as Euler’s rotation theorem. Hence, it is possible to express the rotation between two coordinate frames in terms of an angle $\alpha$ and a unit vector $n$ around which the rotation takes place. In this section, we will derive a relation between the representation $\alpha,n$ and the rotation matrix parametrization from the previous section. Instead of directly considering the rotation of a coordinate frame, we start by considering the rotation of a vector. Note that a counterclockwise rotation of the coordinate frame is equivalent to a clockwise rotation of a vector, see Example 3.2.

Example 3.2** (Rotation of a coordinate frame and rotation of a vector)**

Consider the 2D example in Figure 3.3, where on the left, a vector $x$ is rotated clockwise by an angle $\alpha$ to $x_{\star}$. This is equivalent to (on the right) rotating the coordinate frame $v$ counterclockwise by an angle $\alpha$. Note that $x^{\text{v}}_{\star}=x^{\text{u}}$.

**

In Figure 3.4, a vector $x$ is rotated an angle $\alpha$ around the unit vector $n$. We denote the rotated vector by $x_{\star}$. Suppose that $x$ as expressed in the coordinate frame $v$ is known (and denoted $x^{\text{v}}$) and that we want to express $x_{\star}^{\text{v}}$ in terms of $x^{\text{v}}$, $\alpha$ and $n$. It can first be recognized that the vector $x$ can be decomposed into a component parallel to the axis $n$, denoted $x_{\parallel}$, and a component orthogonal to it, denoted $x_{\perp}$, as

[TABLE]

where $\cdot$ denotes the inner product. Similarly, $x_{\star}^{\text{v}}$ can be decomposed as

[TABLE]

Hence, $x^{\text{v}}_{\star}$ can be expressed in terms of $x^{\text{v}}$ as

[TABLE]

Denoting the rotated coordinate frame the $u$-frame and using the equivalence between $x^{\text{v}}_{\star}$ and $x^{\text{u}}$ as shown in Example 3.2, this implies that

[TABLE]

This equation is commonly referred to as the rotation formula or Euler’s formula [135]. Note that the combination of $n$ and $\alpha$, or $\eta=n\alpha$, is denoted as the rotation vector or the axis-angle parameterization.

To show the equivalence between (3.13) and the rotation matrix parametrization, we will rewrite (3.13). Here, we make use of the fact that a cross product can be written as a matrix vector product. Given vectors $u$ and $v$ we have,

[TABLE]

where $u_{1},u_{2},u_{3}$ denote the three components of the vector $u$. Furthermore, given vectors $u$, $v$ and $w$, multiple cross products can be expanded in terms of the inner product as

[TABLE]

Using these relations, (3.13) can be rewritten as

[TABLE]

Comparing (3.16) and (3.9a), it can be seen that a rotation matrix can be parametrized in terms of $\alpha,n$ as

[TABLE]

Note that equivalently, $R^{\text{uv}}(n^{\text{v}},\alpha)$ can also be written as

[TABLE]

since

[TABLE]

The rotation vector introduced in this section parametrizes the orientation in only three parameters. It is, however, not a unique parametrization since adding $2\pi$ to any angle $\alpha$ results in the same orientation. This is called wrapping. As shown in (3.17) and (3.18), the rotation matrix can straightforwardly be expressed in terms of the axis-angle representation.

3.2.3 Euler angles

Rotation can also be defined as a consecutive rotation around three axes in terms of so-called Euler angles. We use the convention $(z,y,x)$ which first rotates an angle $\psi$ around the $z$-axis, subsequently an angle $\theta$ around the $y$-axis and finally an angle $\phi$ around the $x$-axis. These angles are illustrated in Figure 3.5. Assuming that the $v$-frame is rotated by $(\psi,\theta,\phi)$ with respect to the $u$-frame as illustrated in this figure, the rotation matrix $R^{\text{uv}}$ is given by

[TABLE]

where we make use of the notation introduced in (3.17) and the following definition of the unit vectors

[TABLE]

The $\psi,\theta,\phi$ angles are also often referred to as yaw (or heading), pitch and roll, respectively. Furthermore, roll and pitch together are often referred to as inclination.

Similar to the rotation vector, Euler angles parametrize orientation as a three-dimensional vector. Euler angle representations are not unique descriptions of a rotation for two reasons. First, due to wrapping of the Euler angles, the rotation $(0,0,0)$ is for instance equal to $(0,0,2\pi k)$ for any integer $k$. Furthermore, setting $\theta=\tfrac{\pi}{2}$ in (3.20), leads to

[TABLE]

Hence, only the rotation $\phi-\psi$ can be observed. Because of this, for example the rotations $(\tfrac{\pi}{2},\tfrac{\pi}{2},0)$, $(0,\tfrac{\pi}{2},-\tfrac{\pi}{2})$, $(\pi,\tfrac{\pi}{2},\tfrac{\pi}{2})$ are all three equivalent. This is called gimbal lock [31].

3.2.4 Unit quaternions

A commonly used parametrization of orientation is that of unit quaternions. Quaternions were first introduced by [52] and are widely used in orientation estimation algorithms, see e.g. [76, 59]. A unit quaternion use a $4$-dimensional representation of the orientation according to

[TABLE]

A unit quaternion is not a unique description of an orientation. The reason for this is that if $q$ represents a certain orientation, then $-q$ describes the same orientation.

A rotation can be defined using unit quaternions as

[TABLE]

where $\left(q^{\text{uv}}\right)^{\mathsf{c}}=q^{\text{vu}}$ denotes the quaternion conjugate, defined as

[TABLE]

and $\mathrm{bar}x^{\text{v}}$ denotes the quaternion representation of $x^{\text{v}}$ as

[TABLE]

Note that (3.26) is typically not a unit quaternion. The notation $\odot$ denotes the quaternion multiplication given by

[TABLE]

where

[TABLE]

Using (3.25)–(3.28), (3.24) can be written as

[TABLE]

Comparing (3.2.4) to (3.17), it can be recognized that if we choose

[TABLE]

the two rotation formulations are equivalent since

[TABLE]

Here, we made use of standard trigonometric relations and the fact that since $\|n^{\text{v}}\|_{2}=1$, $n^{\text{v}}\left(n^{\text{v}}\right)^{\mathsf{T}}=\mathcal{I}_{3}+[n^{\text{v}}\times]^{2}$. Hence, it can be concluded that $q^{\text{uv}}$ can be expressed in terms of $\alpha$ and $n^{\text{v}}$ as in (3.30).

Equivalently, $q^{\text{uv}}(n^{\text{v}},\alpha)$ can also be written as

[TABLE]

where

[TABLE]

This leads to

[TABLE]

Note the similarity to (3.18) and (3.2.2). The reason why both rotation matrices and unit quaternions can be described in terms of an exponential of a rotation vector will be discussed in §3.3.1.

3.3 Probabilistic orientation modeling

The four parametrizations of orientation discussed in §3.2 can be used interchangeably. However, the choice of which parametrization to use as states $x_{t}$ in the filtering and smoothing problems introduced in §3.1 has significant impact on the workings of the algorithm. An important reason for this is that estimation algorithms typically assume that the unknown states and parameters are represented in Euclidean space. For instance, they assume that the subtraction of two orientations gives information about the difference in orientation and that the addition of two orientations is again a valid orientation. For the four parametrizations discussed in §3.2, this is generally not true. For instance, due to wrapping and gimbal lock, subtraction of Euler angles and rotation vectors can result in large numbers even in cases when the rotations are similar. Also, addition and subtraction of unit quaternions and rotation matrices do not in general result in a valid rotation. The equality constraints on the norm of unit quaternions and on the determinant and the orthogonality of rotation matrices are typically hard to include in the estimation algorithms. In §3.3.1, we will discuss a method to represent orientation in estimation algorithms that deals with the issues described above. It is frequently used in the algorithms that will be described in Chapter 4. In §3.3.2, we will also discuss some alternative methods to parametrize orientation for estimation purposes.

3.3.1 Linearization

As mentioned in §3.2, the group of rotations in three dimensions is the special orthogonal group $SO(3)$. More specifically, $SO(3)$ is a so-called matrix Lie group. For a discussion on the properties of matrix Lie groups and on the reasons why $SO(3)$ is indeed such a group we refer the reader to e.g. [8]. Since rotations are a matrix Lie group, there exists an exponential map from a corresponding Lie algebra. Using this property, it is possible to represent orientations on $SO(3)$ using unit quaternions or rotation matrices, while orientation deviations are represented using rotation vectors on $\mathbb{R}^{3}$, see e.g. [14]. Hence, we encode an orientation $q^{\text{nb}}_{t}$ in terms of a linearization point parametrized either as a unit quaternion $\tilde{q}^{\text{nb}}_{t}$ or as a rotation matrix $\tilde{R}^{\text{nb}}_{t}$ and an orientation deviation using a rotation vector $\eta_{t}$. Assuming that the orientation deviation is expressed in the body frame $n$,111A similar derivation can be done by assuming an orientation deviation in the navigation frame $b$.

[TABLE]

where analogously to (3.2.4) and (3.2.2),

[TABLE]

For notational convenience, in the remainder we will use the mappings

[TABLE]

which allow us to rewrite (3.35) as

[TABLE]

The reverse mappings are defined as

[TABLE]

where $\log R$ is the standard matrix logarithm. Since we typically assume that $\eta_{t}^{\text{n}}$ is small, we will frequently make use of the following approximations

[TABLE]

The idea briefly outlined in this section is closely related to approaches used to estimate orientation in robotics, see e.g. [47, 48, 14, 8, 38]. It is also related to the so-called multiplicative extended Kalman filter (MEKF) frequently used in aeronautics, see e.g. [93, 28].

3.3.2 Alternative methods

An alternative method to estimate orientation assumes that the states representing the orientation lie on a manifold. This can be done by modeling the orientation and its uncertainty using a spherical distribution which naturally restricts the orientation estimates and their uncertainties to be in $SO(3)$. In recent years, a number of approaches have been proposed to estimate the orientation using these kinds of distributions. For instance, in [77, 44, 45] algorithms are presented to estimate orientation by modeling it using a Bingham distribution.

The difficulties caused by directly using one of the four orientation parametrizations introduced in §3.2 in orientation estimation algorithms is widely recognized. Nevertheless, a large number of approaches directly uses these parametrizations in estimation algorithms. For instance, it is common practice to use unit quaternions in estimation algorithms and to normalize the resulting quaternions each time they loose their normalization, see e.g. [124, 91, 89]. Different approaches to handle the normalization of the quaternions in these algorithms are discussed in [67].

3.4 Measurement models

In the past two sections, we have focused on how orientations can be parametrized. In this section, we will go back to the probabilistic models for the position and orientation estimation problems introduced in §3.1 and provide different measurement models $p(y_{t}\mid x_{t},\theta)$.

3.4.1 Gyroscope measurement models

As discussed in §2.2, the gyroscope measures the angular velocity $\omega_{\text{ib}}^{\text{b}}$ at each time instance $t$. However, as shown in §2.4, its measurements are corrupted by a slowly time-varying bias $\delta_{\omega,t}$ and noise $e_{\omega,t}$. Hence, the gyroscope measurement model is given by

[TABLE]

As was shown in Figure 2.3, the gyroscope measurement noise is quite Gaussian. Because of this, it is typically assumed that $e_{\omega,t}^{\text{b}}\sim\mathcal{N}(0,\Sigma_{\omega})$. If the sensor is properly calibrated, the measurements in the three gyroscope axes are independent. In that case, it can be assumed that

[TABLE]

The gyroscope bias $\delta_{\omega,t}^{\text{b}}$ is slowly time-varying, as discussed in §2.4. There are two conceptually different ways to treat this slowly time-varying bias. One is to treat the bias as a constant parameter, assuming that it typically changes over a longer time period than the time of the experiment. The bias can then either be pre-calibrated in a separate experiment, or it can be considered to be part of the unknown parameters $\theta$ as introduced in §3.1. Alternatively, it can be assumed to be slowly time-varying. This can be justified either by longer experiment times or by shorter bias stability. In the latter case, $\delta_{\omega,t}^{\text{b}}$ can instead be considered as part of the state vector $x_{t}$ and can for instance be modeled as a random walk

[TABLE]

where $e^{\text{b}}_{\delta_{\omega},t}\sim\mathcal{N}(0,\Sigma_{\delta_{\omega},t})$ represents how constant the gyroscope bias actually is.

Modeling the sensor noise and bias is related to the sensor properties. However, there are also modeling choices related to the experiments that can be made. As described in §2.2, the angular velocity $\omega_{\text{ib}}^{\text{b}}$ can be expressed as

[TABLE]

If the sensor does not travel over significant distances as compared to the size of the earth — which is often the case for the applications discussed in Chapter 1 — the navigation frame $n$ can safely be assumed to be stationary. In that case, the transport rate $\omega_{\text{en},t}^{\text{n}}$ is zero. Although the earth rotation $\omega_{\text{ie}}$ as expressed in the body frame $b$ is not constant, its magnitude as compared to the magnitude of the actual measurements is fairly small (see §2.2 and the experimental data presented in Example 2.3). Assuming that the earth rotation is negligible and the navigation frame is stationary leads to the following simplified measurement model

[TABLE]

3.4.2 Accelerometer measurement models

The accelerometer measures the specific force $f^{\text{b}}_{t}$ at each time instance $t$, see also §2.3. As shown in §2.4, the accelerometer measurements are typically assumed to be corrupted by a bias $\delta_{\text{a},t}$ and noise $e_{\text{a},t}$ as

[TABLE]

The accelerometer noise is typically quite Gaussian as was illustrated in Figure 2.4 and can hence be modeled as $e_{\text{a},t}^{\text{b}}\sim\mathcal{N}(0,\Sigma_{\text{a}})$. For a properly calibrated sensor, the covariance matrix $\Sigma_{\text{a}}$ can often be assumed to be diagonal.

The accelerometer bias $\delta_{\text{a},t}^{\text{b}}$ is slowly time-varying. Similar to the gyroscope bias, the accelerometer bias can either be modeled as a constant parameter, or as part of the time-varying state, for instance using a random walk model as in (3.43).

As introduced in §2.3, the specific force measured by the accelerometer is given by

[TABLE]

Assuming that the navigation frame is fixed to the earth frame, we derived a relation for $a_{\text{ii}}^{\text{n}}$ as

[TABLE]

The centrifugal acceleration $\omega_{\text{ie}}^{\text{n}}\times\omega_{\text{ie}}^{\text{n}}\times p^{\text{n}}$ is typically absorbed in the local gravity vector. The magnitude of the Coriolis acceleration is small compared to the magnitude of the accelerometer measurements (see Example 2.2 and the experimental data presented in Example 2.3). Neglecting this term leads to the following simplified measurement model

[TABLE]

Since the accelerometer measures both the local gravity vector and the linear acceleration of the sensor, it provides information both about the change in position and about the inclination of the sensor. For orientation estimation, only the information about the inclination is of concern. Hence, a model for the linear acceleration needs to be made to express the relation between the inclination and the measurements. To model this, it can be recognized that in practice, most accelerometer measurements are dominated by the gravity vector, as illustrated in Example 3.3.

Example 3.3** (Magnitude of a sensor’s linear acceleration)**

Let us consider a 1D example where a sensor has an initial velocity $v_{1}=0~{}\metrepersecond$ and accelerates with $a_{\text{nn}}^{\text{n}}=9.82~{}\metrepersquaresecond$. After $4.51$ seconds, the sensor will have traveled $100$ meters. This is about twice as fast as the world record currently held by Usain Bolt. In fact, humans can reach fairly high accelerations but can only accelerate for a short time. Naturally, cars can accelerate to higher velocities than humans. The sensor in this example has reached a final velocity of $160~{}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{h}$. Even in the case of a car it is therefore unlikely that it can have an acceleration this high for a long period of time.

**

Since the accelerometer measurements are typically dominated by the gravity vector, a commonly used model assumes the linear acceleration to be approximately zero

[TABLE]

Naturally, the model (3.50) is almost never completely true. However, it can often be used as a sufficiently good approximation of reality. Note that the noise term $e_{\text{a},t}^{\text{b}}$ in this case does not only represent the measurement noise, but also the model uncertainty. The model (3.50) can for instance be used in combination with outlier rejection where measurements that clearly violate the assumption that the linear acceleration is zero are disregarded. It is also possible to adapt the noise covariance matrix $\Sigma_{\text{a}}$, depending on the sensor’s acceleration [39, 115]. Furthermore, it is possible to instead model the acceleration based on physical reasoning [85].

3.4.3 Modeling additional information

In this section we will discuss models for the measurements we use to complement the inertial sensor measurements. For orientation estimation we use magnetometers, while for pose estimation we use position measurements.

Magnetometer models

Magnetometers measure the local magnetic field, consisting of both the earth magnetic field and the magnetic field due to the presence of magnetic material. The (local) earth magnetic field is denoted $m^{\text{n}}$ and it is illustrated in Figure LABEL:fig:models-magneticField. Its horizontal component points towards the earth’s magnetic north pole. The ratio between the horizontal and vertical component depends on the location on the earth and can be expressed in terms of the so-called dip angle $\delta$. The dip angle and the magnitude of the earth magnetic field are accurately known from geophysical studies, see e.g. [102].

Assuming that the sensor does not travel over significant distances as compared to the size of the earth, the local earth magnetic field can be modeled as being constant. In case no magnetic material is present in the vicinity of the sensor, orientation information can be deduced from the magnetometer. More specifically, magnetometers are typically used to complement accelerometers to provide information about the sensor heading, i.e. about the orientation around the gravity vector which can not be determined from the accelerometer measurements. Magnetometers provide information about the heading in all locations on the earth except on the magnetic poles, where the local magnetic field $m^{\text{n}}$ is vertical. Orientation can be estimated based on the direction of the magnetic field. The magnitude of the field is irrelevant. Because of this, without loss of generality we model the earth magnetic field as

[TABLE]

i.e. we assume that $\|m^{\text{n}}\|_{2}=1$. Assuming that the magnetometer only measures the local magnetic field, its measurements $y_{\text{m},t}$ can be modeled as

[TABLE]

where $e_{\text{m},t}\sim\mathcal{N}(0,\Sigma_{\text{m}})$. The noise $e_{\text{m},t}$ represents the magnetometer measurement noise as well as the model uncertainty. Note that using the models (3.51) and (3.52), we define the heading with respect to the magnetic north, which is sufficient for most applications. In case we would be interested in navigation with respect to the true north instead, the magnetic declination needs to be taken into account. Since the magnetic declination depends on the location on the earth it is necessary to know the location of the sensor to correct for this.

In practice, the actual magnetic field can differ significantly from the earth magnetic field. In indoor environments, for instance, presence of magnetic material in the structure of buildings and in furniture influences the magnetic field that is measured by the magnetometer. Furthermore, the magnetic field is affected in applications where the magnetometer is mounted in e.g. a vehicle, train or on a robot. In case the magnetic material is rigidly attached to the sensor, the magnetometer can be calibrated for its presence [75, 149, 118, 126]. The presence of magnetic material in the vicinity of the sensor that can not be calibrated for is of major concern for practical applications. Because of this, there is a vast amount of literature on the topic, see e.g. [21, 80, 120].

Note that when using magnetometers for orientation estimation, the presence of magnetic material is typically considered to be an undesired disturbance. However, the presence of magnetic material can also be considered to be a property which can be exploited. This is done in approaches which use the magnetic field as a source of position information [143, 56, 119].

Position information

Position information can be obtained from for instance GNSS or UWB measurements. In this tutorial, we will consider a very basic measurement model where the sensors directly measure the position as

[TABLE]

with $e_{\text{p},t}\sim\mathcal{N}(0,\Sigma_{\text{p}})$.

Many sensors do not measure the position directly. Their measurements can, however, be pre-processed to obtain position estimates and their corresponding covariances [51]. For example, time of arrival measurements can be pre-processed using multilateration techniques. Measurements of this type often contain a significant amount of outliers. The reason is that the signals can be delayed due to multipath or non-line-of-sight (NLOS) conditions. Possible solutions deal with this by doing outlier rejection, by using robust algorithms, see e.g. [158], or by modeling the noise distribution as a non-Gaussian distribution, see e.g. [74, 51, 104]. We will discuss an example of this in more detail in §6.1.

3.5 Choosing the state and modeling its dynamics

The fundamental continuous-time relations that form the basis of our dynamic models are the fact that the position $p^{\text{n}}$, velocity $v^{\text{n}}_{\text{n}}$ and acceleration $a^{\text{n}}_{\text{nn}}$ are related as

[TABLE]

and that the orientation and angular velocity $\omega_{\text{nb},t}^{\text{b}}$ are related as

[TABLE]

depending on orientation parametrization. For a derivation of (3.55), see e.g. [59]. Using an Euler discretization of (3.54) assuming that the acceleration is constant between samples, the dynamics of the position and velocity can be expressed in terms of the acceleration as

[TABLE]

where $T$ is the time between two samples. Similarly, the dynamics of the orientation can be expressed in terms of unit quaternions or rotation matrices as

[TABLE]

Dynamic models describe how the state changes over time. For the problem of position and orientation estimation using inertial sensors, there are two commonly used modeling alternatives for the dynamics [50]. In the first, the state vector $x_{t}$ is chosen to consists of

[TABLE]

The change in position, velocity and orientation states can then be described in terms of the velocity, acceleration and angular velocity states, respectively. The dynamics of the acceleration and the angular velocity can be described in terms of a motion model. Examples of motion models that can be used are a constant acceleration model, which assumes that the dynamics of the acceleration can be described as

[TABLE]

with $w_{\text{a},t}\sim\mathcal{N}(0,\Sigma_{\text{w,a}})$, and a constant angular velocity model, which describes the dynamics of the angular velocity as

[TABLE]

with $w_{\omega,t}\sim\mathcal{N}(0,\Sigma_{\text{w},\omega})$. The process noise terms $w_{\text{a},t}$ and $w_{\omega,t}$ model the assumptions on how constant the acceleration and angular velocity actually are.

Alternatively, the state vector $x_{t}$ can be chosen as

[TABLE]

To describe the dynamics of the states, the inertial measurements can then be used as an input to the dynamic equation (3.5). Hence, the change in position, velocity and orientation is directly modeled in terms of the inertial measurements. In this case, expressions for $a^{\text{n}}_{\text{nn},t}$ and $\omega^{\text{b}}_{\text{nb},t}$ in (3.56) and (3.57) are obtained from the accelerometer measurement model and the gyroscope measurement model, see §3.4. The process noise can explicitly be modeled in terms of the accelerometer measurement noise $e_{\text{a},t}$ and the gyroscope measurement noise $e_{\omega,t}$.

The benefit of using a motion model for the state dynamics is that knowledge about the motion of the sensor can be included in this model. However, it comes at the expense of having a larger state vector. The benefit of using the inertial measurements as an input to the dynamics is that the process noise has the intuitive interpretation of representing the inertial measurement noise. Hence, the latter approach is often used for applications where it is difficult to obtain sensible motion models. Another difference between the two approaches is that changes in the acceleration and angular velocity will have a slightly faster effect on the state when the inertial measurements are used as an input to the dynamics. The reason for this is that the constant acceleration and angular velocity models delay the effects of changes in the dynamics. Because of this, their estimates typically look slightly more smooth.

3.6 Models for the prior

Looking back at §3.1, to solve the smoothing and filtering problems (3.3a) and (3.3b), we have now discussed different models for the dynamics $p(x_{t}\mid x_{t-1},\theta)$ in §3.5 and for the measurements $p(y_{t}\mid x_{t},\theta)$ in §3.4. The remaining distributions to be defined are $p(x_{1}\mid\theta)$ and $p(\theta)$. Note that we typically assume that $x_{1}$ is independent of $\theta$. Hence, in this section we focus on the prior distributions $p(x_{1})$ and $p(\theta)$.

In many cases, there is fairly little prior information available about the parameters $\theta$. It is, however, often possible to indicate reasonable values for the parameters. For example, it is reasonable to assume that the gyroscope bias is fairly small but can be both positive and negative. For instance, for the data presented in Example 2.3, the average bias over the 55-minute data set was $\begin{pmatrix}35.67&54.13&-1.07\end{pmatrix}^{\mathsf{T}}\cdot 10^{-4}~{}\radianpersecond$. If we would assume that in $68\%$ of the cases, the bias is within the bounds $-\sigma_{\delta_{\omega}}\leq\delta_{\omega,t}^{\text{b}}\leq\sigma_{\delta_{\omega}}$ with $\sigma_{\delta_{\omega}}=5\cdot 10^{-3}$, a reasonable prior would be

[TABLE]

For the prior $p(x_{1})$, it is typically possible to get a reasonable estimate from data. For the position and velocity, this can be modeled as

[TABLE]

Here, the estimate $\breve{p}_{1}^{\text{n}}$ can for instance be determined based on the first position measurement. In that case, the uncertainty $\sigma_{\text{p,i}}$ can also be chosen equal to the uncertainty of the position measurements. In case no additional information is available, the estimates $\breve{p}_{1}^{\text{n}}$ and $\breve{v}_{1}^{\text{n}}$ can be set to zero with an appropriate standard deviation instead.

A commonly used method to determine the initial orientation is to use the first accelerometer and magnetometer samples. This method is based on the fact that given two (or more) linearly independent vectors in two coordinate frames, the rotation between the two coordinate frames can be determined. The implicit assumption is that the accelerometer only measures the gravity vector and the magnetometer only measures the local magnetic field. Hence, the four vectors are given by the measurements $y_{\text{a},t}$ and $y_{\text{m},t}$, the local gravity vector $g^{\text{n}}$ and the local magnetic field $m^{\text{n}}$. These vectors are linearly independent except when the measurements are obtained on the magnetic north or south poles where the dip angle is $\delta=0$ and the magnetic field does not contain any heading information.

The accelerometer provides information about the sensor’s inclination. Heading information is provided by the magnetometer. However, at all locations except on the equator, the magnetometer also provides information about the inclination due to its non-zero vertical component, see (3.51). In practice, the accelerometer typically provides more accurate inclination information. Hence, we choose to use the magnetometer only to provide heading information by projecting the magnetic field and the magnetometer measurement on the horizontal plane. Furthermore, we normalize the vectors. Because of this, an adapted model uses the four normalized vectors

[TABLE]

A number of algorithms are available to estimate the orientation from these vectors. Well-known examples are the TRIAD algorithm, the QUEST algorithm, see e.g. [137], and the method presented in [62]. For our problem at hand, these methods give equivalent results, even though they use slightly different solution strategies. Generally speaking, they solve the problem of determining the rotation $q^{\text{nb}}$ from

[TABLE]

Recall from (3.24) that $q^{\text{nb}}\odot\mathrm{bar}{x}^{\text{b}}\odot q^{\text{bn}}$ is the rotation of the vector $x^{\text{b}}$ to the $n$-frame. The optimization problem (3.6) therefore determines the orientation $q^{\text{nb}}$ that minimizes the distance between the normalized magnetic field and gravity vectors measured in the first sample and the normalized magnetic field and gravity vectors in the navigation frame. These four vectors were defined in (3.64).

Defining

[TABLE]

where the left and right quaternion multiplications are defined in (3.28), (3.6) can equivalently be written as

[TABLE]

For a derivation, see [59]. The solution to this problem is given by the eigenvector corresponding to the largest eigenvalue of $A$. Note that although this method can be used to compute the orientation from any two linearly independent vectors in two coordinate frames, we only use it to compute a prior on the orientation.

Based on the estimate $\breve{q}_{1}^{\text{nb}}$ from (3.67), we can model the orientation at time $t=1$ in terms of an orientation deviation

[TABLE]

Explicit formulations for the covariance of the orientation estimates from the TRIAD and QUEST algorithms are discussed by [136]. In practice, however, the accuracy of the estimates from (3.67) highly depends on the validity of the model assumptions, i.e. on whether the sensor is indeed far from magnetic material and whether the linear acceleration is indeed zero. Because this has such a significant influence on the quality of the estimates, we choose $\Sigma_{\eta,\text{i}}$ and $\Sigma_{\text{q,i}}$ somewhat conservatively. Modeling that in $68\%$ of the cases the orientation error is less than $20^{\circ}$,

[TABLE]

where we use the fact that $(q^{\mathsf{R}})^{\mathsf{T}}=(q^{\mathsf{c}})^{\mathsf{R}}$. Note that the covariance $\Sigma_{\text{q,i}}$ is defined in terms of the covariance $\Sigma_{\eta,\text{i}}$ to allow for explicit comparison between different algorithms in Chapter 4.

3.7 Resulting probabilistic models

The information from the previous sections can now be combined into one model which will be used in the algorithms in Chapter 4. In this section, we describe our modeling choices for the pose estimation problem and for the orientation estimation problem.

We assume that the sensor does not travel over significant distances as compared to the size of the earth and hence keep the navigation frame $n$ fixed with respect to the earth frame $e$. Furthermore, we assume that the magnitude of the earth rotation and of the Coriolis acceleration are negligible. Our gyroscope and accelerometer models are hence given by

[TABLE]

In the remainder, for notational convenience we drop the subscripts $n$ which indicate in which frame the differentiation is performed, see §2.3, and use the shorthand notation $a^{\text{n}}$ for $a_{\text{nn}}^{\text{n}}$. Furthermore, we will denote $\omega_{\text{nb}}^{\text{b}}$ simply by $\omega$ and omit the superscript $b$ on the noise terms $e_{\text{a},t}$ and $e_{\omega,t}$ and the bias terms $\delta_{\text{a},t}$ and $\delta_{\omega,t}$. We assume that the inertial measurement noise is given by

[TABLE]

i.e. we assume that the three sensor axes are independent and have the same noise levels.

3.7.1 Pose estimation

For pose estimation, we model the accelerometer and gyroscope measurements as inputs to the dynamics. Hence, the state vector consists of the position $p^{\text{n}}_{t}$, the velocity $v^{\text{n}}_{t}$ and a parametrization of the orientation. We use the inertial measurements in combination with position measurements to estimate the pose.

Using the accelerometer measurement model (3.70a) in (3.56), the dynamics of the position and velocity is given by

[TABLE]

where without loss of generality, we switch the sign on the noise. Note that the noise term $e_{\text{a},t}$ should be rotated to the navigation frame $n$ by multiplying it with the rotation matrix $R^{\text{nb}}_{t}$. However, because of the assumption (3.71), the rotation matrix can be omitted without loss of generality. The dynamics of the orientation parametrized using quaternions is given by

[TABLE]

Equivalent dynamic models can be obtained for the other parametrizations.

The position measurements are modeled as in (3.53). In summary, this leads to the following state space model for pose estimation

[TABLE]

where the initial orientation uncertainty is given in terms of a standard deviation of $20^{\circ}$.

In Chapter 4 we assume that the inertial measurement are properly calibrated. Hence, we assume that their biases $\delta_{\text{a},t}$ and $\delta_{\omega,t}$ are zero. Calibration is the topic of Chapter 5 where we will also introduce possible extensions of the state space model in which the bias terms are included either as states or as unknown parameters.

3.7.2 Orientation estimation

For orientation estimation, the state vector only consists of a parametrization of the orientation. We use the inertial sensors in combination with the magnetometer measurements to estimate the orientation. The magnetometer measurements are modeled as in (3.52). Instead of using the accelerometer measurements model (3.70a), we use the model (3.50) where it is assumed that the linear acceleration is approximately zero. This leads to the following state space model for orientation estimation,

[TABLE]

with $\Sigma_{\omega}=\sigma_{\omega}^{2}\,\mathcal{I}_{3}$ and $\Sigma_{\text{a}}=\sigma_{\text{a}}^{2}\,\mathcal{I}_{3}$. The initial orientation is given by the QUEST algorithm described in §3.6 and is modeled as in (3.74g) or (3.74h). Also for orientation estimation, in Chapter 4 we assume that the inertial measurements are properly calibrated. Hence, we assume that the bias $\delta_{\omega,t}$ is zero.

Chapter 4 Estimating Position and Orientation

In this chapter we will focus on position and orientation estimation using the models (3.74) and (3.75) derived in Chapter 3. In §4.1, we will first describe a method to solve the smoothing problem (3.3a). Subsequently, in §4.2–§4.4, we will derive different methods for solving the filtering problem (3.3b). In each section, after a general introduction of the estimation method, we will illustrate the method by explicitly deriving algorithms to estimate the orientation using the state space model (3.75). The orientation estimation problem also illustrates the most important parts of the pose estimation problem, since most complexities lie in the parametrization of the orientation and in the nonlinear nature of the orientation. In §4.5, we show some characteristics of the different algorithms for the orientation estimation problem. In §4.6, we will discuss how the algorithms for orientation estimation can be extended to also estimate the position. Throughout this section, we assume that the sensors are calibrated, i.e. we assume that we do not have any unknown parameters $\theta$ in our models. Because of this, the models that we use are the most basic models that can be used for position and orientation estimation using inertial sensors.

4.1 Smoothing in an optimization framework

Perhaps the most intuitive way to solve the smoothing problem is by posing it as an optimization problem, where a maximum a posteriori (MAP) estimate is obtained as

[TABLE]

Here, we use the notation in terms of probability distributions as introduced in §3.1 and model the measurements and the state dynamics as described in §3.4 and §3.5, respectively. Furthermore, we assume that a prior on the initial state is obtained using the measurements at $t=1$ as described in §3.6. Because of this, the measurement model $p(y_{1}\mid x_{1})$ from (3.3a) is explicitly omitted in (4.1). Note that in practice, we typically minimize $-\log p(x_{1:N}\mid y_{1:N})$ instead of maximizing $p(x_{1:N}\mid y_{1:N})$ itself, resulting in the optimization problem

[TABLE]

There are various ways to solve problems of this kind, for instance particle smoothers [81], an extended Rauch-Tung-Striebel (RTS) smoother [127] and optimization methods, see e.g. [103, 95]. The latter approach is closely related to iterated Kalman smoothers [10, 65]. We will solve the problem using an optimization method. Compared to extended RTS smoothers, optimization methods allow for more flexibility in the models that are being used. For instance, additional information outside of the standard state space model can straightforwardly be included. Optimization approaches are typically computationally heavier than extended RTS smoothers but less heavy than particle smoothers. The latter are capable of capturing the whole distribution, which is a clear advantage when the distributions are multi-modal. Optimization instead gives a point estimate and an associated measure of uncertainty. This is typically sufficient for position and orientation estimation using inertial sensors.

4.1.1 Gauss-Newton optimization

To obtain a smoothing estimate of the position and orientation using optimization, we first recognize that for our models (3.74) and (3.75), all probability distributions in (4.2) are Gaussian. Let us therefore consider a slightly more general problem where the objective function consists of the product of Gaussian probability functions $p(e_{i}(x_{1:N}))$, $i=1,\ldots,M$. Hence, the optimization problem can be written as

[TABLE]

The probability distribution of $e_{i}(x)$ is given by

[TABLE]

Omitting the terms independent of $x_{1:N}$, the optimization problem (4.3) reduces to

[TABLE]

with $\|e_{i}(x_{1:N})\|_{\Sigma_{i}^{-1}}^{2}=e^{\mathsf{T}}_{i}(x_{1:N})\Sigma_{i}^{-1}e_{i}(x_{1:N})$. The function that is being minimized in optimization problems, is often referred to as the objective function.

The solution to (4.5) can be found by studying the shape of the objective function as a function of $x_{1:N}$. This can be characterized in terms of the gradient $\mathcal{G}(x_{1:N})$ and Hessian $\mathcal{H}(x_{1:N})$, which provide information about the slope and curvature of the function, respectively. Defining

[TABLE]

and the stacked variables

[TABLE]

the gradient and the Hessian are given by

[TABLE]

Note that for notational convenience, we have omitted the explicit dependence of $\varepsilon$ on $x_{1:N}$. In (4.6), we introduced the notation $\mathcal{J}(x_{1:N})$, which is the Jacobian of the vector $\varepsilon$ with respect to $x_{1:N}$ as

[TABLE]

where $n_{\varepsilon}$ is the length of the vector $\varepsilon_{i}$. Instead of computing the true Hessian (4.6b), we compute an approximation of it [103], given by

[TABLE]

This has the benefit of not having to compute second derivatives, at the same time as it guarantees that the Hessian is positive semidefinite. The downside of using (4.8) is that it introduces an approximation.

The gradient and the (approximate) Hessian can be used to find the minimum of the objective function. For our models (3.74) and (3.75), in which the functions $e_{i}(x_{1:N})$ are nonlinear, an estimate $\hat{x}_{1:N}$ can iteratively be computed as

[TABLE]

where $k$ denotes the iteration number. The step length $\beta^{(k)}$ is computed for instance using a backtracking line search [103, 16]. The search direction is computed as $\left(\mathcal{\hat{H}}(\hat{x}_{1:N}^{(k)})\right)^{-1}\mathcal{G}(\hat{x}_{1:N}^{(k)})$. Note that an initial point $\hat{x}_{1:N}^{(0)}$ needs to be chosen close enough to the desired minimum to ensure convergence to this minimum.

In case the functions $e_{i}(x_{1:N})$ would be linear, the problem (4.5) would be a least squares (LS) problem for which the update (4.9) directly leads to the minimum of the objective function, irrespective of the initial point. In our case where the functions $e_{i}(x_{1:N})$ are nonlinear due to the nonlinear nature of the orientation, the problem (4.5) is instead a nonlinear least squares (NLS) problem. Each iteration in (4.9) can be interpreted as solving a LS problem around a linearization point. The linearization point is updated after each iteration, bringing it closer to the minimum of the objective function. Computing an estimate $\hat{x}_{1:N}$ by iterating (4.9) until convergence, making use of the approximate Hessian from (4.8), is called Gauss-Newton optimization.

Note that since the inertial sensors sample at high sampling rates, the length of the vector $x_{1:N}$ quickly becomes fairly large. For instance, for inertial sensors sampling at $100\mathrm{Hz}$ for $10$ seconds, $N=1\thinspace 000$ and the size of the (approximate) Hessian $\mathcal{\hat{H}}(x)$ is $1\thinspace 000n_{x}\times 1\thinspace 000n_{x}$, where $n_{x}$ is the length of the vector $x_{t}$. However, as can be seen in (4.2), the components of the objective function only depend on the current and next time steps $x_{t}$ and $x_{t+1}$. Hence, the structure of the Hessian (4.8) is of the form given in Figure 4.1. There exist efficient algorithms to compute search directions for problems with this sparsity pattern, which can be exploited using sparse matrix packages, see e.g. [29], or by using tools like dynamic programming and message passing [11, 46, 123].

4.1.2 Smoothing estimates of the orientation using optimization

In this section, we will illustrate the use of Gauss-Newton optimization to obtain smoothing estimates of the orientation. As discussed in the previous section, the crucial part is to identify the objective function and its Jacobian. From this, the gradient and approximate Hessian can be computed using (4.6a) and (4.8), which can be used to iteratively update the estimates using (4.9).

Combining the general formulation of the smoothing problem (4.2) and using the model for orientation estimation (3.75), the orientation smoothing problem is given by

[TABLE]

with

[TABLE]

Convex equality constraints can straightforwardly be incorporated in optimization problems. However, using $N$ equality constraints to preserve the unit norm of the quaternions, severely complicates the optimization problem since these norm constraints are non-convex. Instead, we encode an orientation in terms of a linearization point parametrized as a unit quaternion $\tilde{q}^{\text{nb}}_{t}$ and an orientation deviation parametrized as a rotation vector $\eta_{t}^{\text{n}}$ as discussed in §3.3.1. Hence, we model the orientation as

[TABLE]

At each Gauss-Newton iteration (4.9), we estimate the state vector $\eta_{1:N}^{\text{n}}$. Before starting the next iteration, the linearization points $\tilde{q}^{\text{nb}}_{1:N}$ are updated and the state vector $\eta_{1:N}^{\text{n}}$ is reset to zero.

Using the notation introduced in §3.3.1, the objective function (4.10) can be expressed in terms of the orientation deviation $\eta_{t}^{\text{n}}$ by rewriting (4.11) as

[TABLE]

Here, $\tilde{R}^{\text{bn}}_{t}$ is the rotation matrix representation of the linearization point $\tilde{q}^{\text{bn}}_{t}$. This leads to the following derivatives

[TABLE]

where, using (3.40) and the definition of the quaternion conjugate (3.25),

[TABLE]

Using the approximate derivatives (4.14), the gradient and approximate Hessian can be computed for the Gauss-Newton iterations. The resulting solution is summarized in Algorithm 1. One of the inputs to the algorithm is an initial estimate of the orientation $\tilde{q}^{\text{nb},(0)}_{1:N}$, which will aid the convergence to the desired minimum. There are (at least) two ways to obtain good initial estimates $\tilde{q}^{\text{nb},(0)}_{1:N}$. First, they can be obtained by direct integration of the gyroscope measurements. As discussed in §1.2, these estimates will suffer from integration drift. However, they still form a good initial estimate for the optimization problem. Second, one of the other, computationally cheaper estimation algorithms that will be discussed in the remainder of this section can be used to obtain initial estimates of the orientation.

4.1.3 Computing the uncertainty

As discussed in §3.1, we are not only interested in obtaining point estimates of the position and orientation, but also in estimating their uncertainty. In our case of Gaussian noise, this is characterized by the covariance. As shown in e.g. [82, 150], if our position and orientation estimation problems would be LS problems, the covariance of the estimates would be given by the inverse of the Hessian of the objective function (4.6b). Instead, we solve an NLS problem, for which a number of LS problems are solved around linearization points closer and closer to the minimum of the objective function. Hence, when the algorithm has converged, the problem can locally well be described by the quadratic approximation around its resulting linearization point. We can therefore approximate the covariance of our estimates as

[TABLE]

An intuition behind this expression is that the accuracy of the estimates is related to the sensitivity of the objective function with respect to the states. The covariance of the estimate will play a crucial role in the filtering approaches discussed in §4.2 and §4.3.

The matrix $\mathcal{J}^{\mathsf{T}}(\hat{x}_{1:N})\mathcal{J}(\hat{x}_{1:N})$ quickly becomes fairly large due to the high sampling rates of the inertial sensors. Hence, computing its inverse can be computationally costly. We are, however, typically only interested in a subset of the inverse. For instance, we are often only interested in diagonal or block diagonal elements representing $\operatorname{cov}(x_{t})$. It is therefore not necessary to explicitly form the complete inverse, which makes the computation tractable also for larger problem sizes.

4.2 Filtering estimation in an optimization framework

One of the benefits of using a smoothing formulation is that all measurements $y_{1:N}$ are used to get the best estimates of the states $x_{1:N}$. However, both the computational cost and the memory requirements grow with the length of the data set. Furthermore, it is a post-processing solution in which we have to wait until all data is available. Alternatively, the filtering problem can be formulated as

[TABLE]

Note that without loss of generality, we have shifted our time indices as compared to the notation in §3.1. The probability distribution $p(x_{t+1}\mid y_{1:t})$ can now be seen as a prior and can be obtained by marginalizing out the previous state $x_{t}$ as

[TABLE]

Assuming that

[TABLE]

the integral in (4.2) can be approximated as

[TABLE]

Here, F_{t}=\left.{\tfrac{\partial f(x_{t})}{\partial x_{t}}}_{\leavevmode\hbox{\set@color\leavevmode\hbox{\set@color}\raisebox{-1.0pt}{\leavevmode\hbox{\set@color}}}}\right|_{\leavevmode\hbox{\set@color\leavevmode\hbox{\set@color\scriptscriptstyle x_{t}=\hat{x}{t}}\raisebox{9.08334pt}{\leavevmode\hbox{\set@color\scriptscriptstyle w{t}=0}}\hskip 12.30035pt}}, \left.{G_{t}=\tfrac{\partial f(x_{t})}{\partial w_{t}}}_{\leavevmode\hbox{\set@color\leavevmode\hbox{\set@color}\raisebox{-1.0pt}{\leavevmode\hbox{\set@color}}}}\right|_{\leavevmode\hbox{\set@color\leavevmode\hbox{\set@color\scriptscriptstyle x_{t}=\hat{x}{t}}\raisebox{9.08334pt}{\leavevmode\hbox{\set@color\scriptscriptstyle w{t}=0}}\hskip 12.30035pt}} and $w_{t}$ is the process noise defined in (3.5). Using (4.22), the filtering problem (4.2) can again be written as a nonlinear least squares problem and can be solved using Gauss-Newton optimization as described in §4.1.1. Instead of optimizing over the whole data set at once as in §4.1, we solve a small optimization problem at each time instance $t$. This is closely related to what is often called an iterated extended Kalman filter [65, 141].

4.2.1 Filtering estimates of the orientation using optimization

In this section, we will derive an algorithm to obtain filtering estimates of the orientation using Gauss-Newton optimization. The resulting algorithm is closely related to the iterated multiplicative extended Kalman filter presented in [142].

To derive the approximation (4.22) for the orientation estimation problem, we first derive an expression for the function $f$ that expresses $\eta_{t+1}^{\text{n}}$ in terms of $\eta^{\text{n}}_{t}$. Using (4.12) in combination with (3.75a),

[TABLE]

Here, $\eta^{\text{n}}_{t}$ and $\tilde{q}_{t}^{\text{nb}}$ provide information about the previous time step. We denote the orientation estimate from the previous time step as $\hat{q}_{t}^{\text{nb}}$ and assume that $\hat{\eta}_{t}^{\text{n}}$ is zero since the linearization point is updated after each iteration of the optimization algorithm. To start with a fairly good linearization point, we choose $\tilde{q}_{t+1}^{\text{nb},(0)}$ at iteration $k=0$ as

[TABLE]

Following (4.22), the distribution $p(\eta^{\text{n}}_{t+1}\mid y_{1:t})$ around this linearization point can now be written as

[TABLE]

with

[TABLE]

The covariance $P_{t\mid t}$ in (4.25) can be approximated as the inverse of the Hessian of the objective function from the previous time step, see also §4.1.3. We make use of the shorthand notation $P_{t+1\mid t}=F_{t}P_{t\mid t}F_{t}^{\mathsf{T}}+G_{t}QG_{t}^{\mathsf{T}}$ and define

[TABLE]

Note that $e_{\text{f},t}$ is equal to zero for iteration $k=0$ but can be non-zero for subsequent iterations. Using this notation, the filtering problem (4.2) results in the following optimization problem

[TABLE]

Note the similarity of this optimization problem to the smoothing formulation in (4.10). The term $\|e_{\text{f},t}\|_{P_{t\mid t-1}^{-1}}$ takes into account both the knowledge about the previous state $x_{t}$ and the dynamics. Furthermore, due to the fact that $P_{t\mid t-1}^{-1}$ is time-varying, the uncertainty and cross-correlation of the states at the previous time instance is taken into consideration. Including this term is similar to the inclusion of an arrival cost in moving horizon estimation approaches [114].

After each Gauss-Newton iteration, we need to recompute the linearization point as

[TABLE]

We also need to compute the covariance around this updated linearization point as

[TABLE]

where $J_{t}^{(k)}$ can be derived similarly to the derivation of $F_{t}$ in (4.26a). Since $J_{t}^{(k)}=\mathcal{I}_{3}$, in practice the relinearization of the covariance can be omitted. The process of estimating orientation using filtering is summarized in Algorithm 2.

4.3 Extended Kalman filtering

Instead of using optimization for position and orientation estimation, it is also possible to use extended Kalman filtering. Extended Kalman filters (EKFs) compute filtering estimates in terms of the conditional probability distribution (3.2). Hence, the approach is similar to the one discussed in §4.2. In fact, extended Kalman filtering can be interpreted as Gauss-Newton optimization of the filtering problem using only one iteration (4.9) with a step length of one, see e.g. [141]. In this section we first introduce how an EKF can be used to compute state estimates in a general nonlinear state space model. Subsequently, we illustrate the use of EKFs for position and orientation estimation by focusing on the orientation estimation problem. Two different implementations will be discussed. First, we introduce an EKF implementation that uses unit quaternions as states. Subsequently, we discuss an EKF which parametrizes the orientation in terms of an orientation deviation from a linearization point, similar to the approach used in §4.1 and §4.2.

An EKF makes use of a nonlinear state space model. Assuming that the measurement noise is additive and that both the process and the measurement noise are zero-mean Gaussian with constant covariance, the state space model is given by

[TABLE]

with $w_{t}\sim\mathcal{N}(0,Q)$ and $e_{t}\sim\mathcal{N}(0,R)$.

The state is estimated recursively by performing a time update and a measurement update. The time update uses the model (4.34a) to “predict” the state to the next time step according to

[TABLE]

with

[TABLE]

Here, the matrix $P$ denotes the state covariance. The double subscripts on $\hat{x}_{t+1\mid t}$ and $P_{t+1\mid t}$ denote the state estimate and the state covariance at time $t+1$ given measurements up to time $t$. Similarly, $\hat{x}_{t\mid t}$ and $P_{t\mid t}$ denote the state estimate and the state covariance at time $t$ given measurements up to time $t$.

The measurement update makes use of the measurement model (4.34b) in combination with the measurements $y_{t}$ to update the “predicted” state estimate as

[TABLE]

with

[TABLE]

and

[TABLE]

Note that in (4.37) we have shifted our notation one time step compared to the notation in (4.35) to avoid cluttering the notation. The EKF iteratively performs a time update and a measurement update to estimate the state and its covariance.

4.3.1 Estimating orientation using quaternions as states

In this section, we will illustrate the use of an EKF to compute filtering estimates of the orientation. As discussed in the previous section, the crucial part is to compute the matrices $F_{t}$, $G_{t}$ and $H_{t}$ to perform the EKF time and measurement updates. Using the state space model (3.75) and using unit quaternions as states in the EKF, the dynamic model is given by

[TABLE]

The following derivatives of the dynamic model (4.40) can be obtained

[TABLE]

In the measurement update of the EKF, the state is updated using the accelerometer and magnetometer measurements. Using the measurement models

[TABLE]

the matrix $H_{t}$ in (4.39) is given by

[TABLE]

The derivative $\tfrac{\partial R^{\text{bn}}_{t\mid t-1}}{\partial q^{\text{nb}}_{t\mid t-1}}$ can be computed from the definition of the relation between and the rotation matrix and the quaternion representation given in (A.10).

Note that the quaternion obtained from the measurement update (4.37) is no longer normalized. We denote this unnormalized quaternion and its covariance by $\tilde{q}_{t\mid t}^{\text{nb}}$ and $\tilde{P}_{t\mid t}$, respectively. The $\tilde{\cdot}$ instead of the $\hat{\cdot}$ is meant to explicitly indicate that the quaternion still needs to be updated. This is done by an additional renormalization step as compared to a standard EKF implementation. A possible interpretation of the renormalization is as an additional measurement update without measurement noise. Hence, we adjust both the quaternion and its covariance as

[TABLE]

Here, $q^{\mathsf{T}}$ is the standard vector transpose. The resulting EKF is summarized in Algorithm 3.

4.3.2 Estimating orientation using orientation deviations as states

An alternative EKF implementation parametrizes the orientation in terms of an orientation deviation around a linearization point. The linearization point is parametrized in terms of quaternions or rotation matrices and denoted $\tilde{q}^{\text{nb}}_{t}$ or equivalently $\tilde{R}^{\text{nb}}_{t}$. The orientation deviation $\eta_{t}^{\text{n}}$ is the state vector in the EKF. This EKF implementation is sometimes referred to as a multiplicative EKF [28, 93]. One of its advantages is that its implementation is computationally attractive since it only uses a 3-dimensional state compared to the 4-dimensional state in Algorithm 3.

Similarly to (4.2.1) in §4.2, the dynamics of the state $\eta^{\text{n}}$ is given by

[TABLE]

In the time update of the EKF, we first use (3.75a) to directly update the linearization point as

[TABLE]

Around this linearization point, the derivatives of the dynamic model (4.48) are given by

[TABLE]

Note the similarity with (4.26).

In the measurement update of the EKF, the state $\eta^{\text{n}}_{t}$ is updated using the accelerometer and magnetometer measurements. The accelerometer measurement equation can be written in terms of the state $\eta^{\text{n}}_{t}$ as

[TABLE]

Equivalently, the magnetometer measurement equation can be written in terms of the state $\eta_{t}^{\text{n}}$ as

[TABLE]

From these equations, the derivatives $H_{t}$ as defined in (4.39) can straightforwardly be computed.

After the measurement update, the orientation deviation $\hat{\eta}^{\text{n}}_{t}$ is non-zero. Hence, as an additional step in the EKF, we update the linearization point and reset the state. In our algorithm, we consider the relinearization as the “measurement update” for the linearization point, i.e. the relinearization updates the linearization point $\tilde{q}^{\text{nb}}_{t\mid t-1}$ to $\tilde{q}^{\text{nb}}_{t\mid t}$ as

[TABLE]

For the same reasons as in §4.2, the relinearization of the covariance can be omitted. The resulting EKF is summarized in Algorithm 4. Note the similarities between this algorithm and Algorithm 2.

4.4 Complementary filtering

An alternative to using EKFs for orientation estimation is to use complementary filters [58, 17]. This type of algorithms again estimates the orientation at time $t$ given the measurements $y_{1:t}$. However, it does not use the probabilistic models presented in Chapter 3. Instead, complementary filters explicitly use the fact that both the gyroscope and the combination of the accelerometer and the magnetometer provide information about the orientation of the sensor. The orientation estimates obtained from the accelerometer and the magnetometer measurements are noisy but accurate over long periods of time. On the other hand, the orientation estimates using the gyroscope measurements are accurate on a short time scale but they drift over time. These properties can be interpreted and exploited in the frequency domain. The orientation estimates using the gyroscope have desirable properties at high frequencies. We would therefore like to filter these using a high-pass filter. Conversely, the orientation estimates obtained from the accelerometer and magnetometer measurements have desirable properties at low frequencies. We would therefore like to filter these orientation estimates using a low-pass filter.

This can be illustrated by considering the one-dimensional case of estimating an angle $\theta$ from gyroscope measurements $y_{\omega}$ and magnetometer measurements $y_{\text{m}}$. We denote the angle estimated by the complementary filter by $\hat{\theta}$, the angle obtained from the magnetometer measurements by $\theta_{\text{m}}$ and the angle obtained from the gyroscope measurements by $\theta_{\omega}$. Note that the latter is obtained by integrating the gyroscope measurements. The Laplace transforms of $\theta$, $\theta_{\text{m}}$, $\theta_{\omega}$ and $y_{\omega}$ are denoted by $\Theta(s)$, $\Theta_{\text{m}}(s)$, $\Theta_{\omega}(s)$ and $Y_{\omega}(s)$, respectively. The complementary filter computes $\Theta(s)$ as

[TABLE]

where $G(s)$ is a low-pass filter and $1-G(s)$ is hence a high-pass filter. Note that the sum of the two filters, $G(s)$ and $1-G(s)$, is equal to one, which is the reason behind the name complementary filter. Choosing $G(s)$ as a first-order low-pass filter $G(s)=\tfrac{1}{as+1}$ and using Euler backward discretization, the filter (4.4) can be written in discrete time as

[TABLE]

where $\gamma=\tfrac{a}{T+a}$. The relation (4.58) allows us to recursively estimate the angle $\hat{\theta}$. The only parameter that needs to be chosen to run the complementary filter (4.58) is the parameter $a$ of the low-pass filter $G(s)$. Choosing a large $a$ (and hence a $\gamma$ close to one) results in a lower cut-off frequency and a more significant contribution of the gyroscope measurements on the orientation estimates. Conversely, a small $a$ (and hence a $\gamma$ close to zero) results in a higher cut-off frequency and a more significant contribution of the magnetometer measurements on the orientation estimates.

There is a strong relationship between complementary and Kalman filtering for linear models, see e.g. [58, 17]. To highlight the similarities and differences between complementary and extended Kalman filtering for orientation estimation, let us define the estimate from the complementary filter as $\hat{\theta}_{t\mid t}$ and write the recursion (4.58) equivalently as

[TABLE]

So far, we have discussed the one-dimensional case of estimating an angle $\theta$. Let us now focus on estimating orientation in three dimensions and parametrize this as a unit quaternion $q^{\text{nb}}$. Two well-known complementary filter implementations to estimate the orientation $q^{\text{nb}}$ are presented in [90] and [89]. Open-source implementations of both algorithms are available online [155]. The filter presented in [89] is specifically designed to be computationally efficient and has for instance been implemented in the robot operating system (ROS) [146]. In this section, we will derive a complementary filter that is inspired by this implementation but which is also meant to illustrate the similarities and differences between complementary filters and extended Kalman filters for orientation estimation. Because of that, the complementary filter that we will present is not as computationally efficient as the one in [89].

Similar to the one-dimensional case (4.59a), the orientation estimate of the complementary filter can be updated using the gyroscope measurements as

[TABLE]

where $\hat{q}^{\text{nb}}$ is the orientation estimate from the complementary filter and the double subscripts are defined analogously to (4.59). Note that the update (4.60) is equivalent to the time update in the EKF in Algorithm 3.

The orientation from the accelerometer and magnetometer measurements, denoted $R_{\text{am},t}^{\text{nb}}$, can be obtained by solving

[TABLE]

see also §3.6. Similar to [89], we do not completely solve the optimization problem (4.4) in each recursion of the complementary filter. Instead, we use an initial estimate $R_{\text{am},t}^{\text{nb},(0)}=\hat{q}^{\text{nb}}_{t\mid t-1}$ and only perform one iteration of an optimization algorithm. In [89], one iteration of a gradient descent algorithm is performed. We instead perform one Gauss-Newton update as in (4.9). Hence, using

[TABLE]

the estimate $q_{\text{am},t}^{\text{nb}}$ is given by

[TABLE]

Analogously to (4.59b), the second update of the complementary filter is therefore given by

[TABLE]

There are many similarities between (4.4) and the measurement update in the EKF in Algorithm 3. First of all, if we neglect the presence of $\Sigma_{\text{a}}$ and $\Sigma_{\text{m}}$ for a moment, $\mathcal{J}_{t}$ and $\varepsilon_{t}$ in (4.62) are equal to $-H_{t}$ and $\varepsilon_{t}$ in the measurement update in Algorithm 3. Inclusion of $\Sigma_{\text{a}}$ and $\Sigma_{\text{m}}$ in (4.62) ensures that the uncertainty of the accelerometer and magnetometer measurements is properly taken into account when obtaining $q_{\text{am},t}^{\text{nb}}$. The contributions to the orientation estimate $\hat{q}_{t\mid t}^{\text{nb}}$ from integration of the gyroscope measurements and from $q_{\text{am},t}^{\text{nb}}$ are determined by the factor $(1-\gamma)\beta$. In the EKF on the other hand, the uncertainty of the gyroscope, accelerometer and magnetometer measurements is taken into account through the covariance matrices $\Sigma_{\omega}$, $\Sigma_{\text{a}}$ and $\Sigma_{\text{m}}$.

Second of all, comparing (4.4) to the EKF measurement update (4.37), the update (4.4) can be interpreted as a scaled version of an EKF update with $P_{t\mid t-1}=\mathcal{I}_{4}$ and $R=0_{3\times 3}$. In that case, the matrix $K_{t}$ in the EKF would be given by $K_{t}=H_{t}^{\mathsf{T}}\left(H_{t}H_{t}^{\mathsf{T}}\right)^{-1}$. The matrix $H_{t}H_{t}^{\mathsf{T}}$ is rank deficient. Because of this, to compute $K_{t}$, we would have to use the pseudo-inverse, see e.g. [46]. Denoting the pseudo-inverse of a matrix as $\cdot^{\dagger}$, $K_{t}=H_{t}^{\mathsf{T}}\left(H_{t}H_{t}^{\mathsf{T}}\right)^{\dagger}=\left(H_{t}^{\mathsf{T}}H_{t}\right)^{-1}H_{t}^{\mathsf{T}}$. An important difference between the EKF and complementary update is that the scaling factor $(1-\gamma)\beta$ in (4.4) is constant, while the matrix $P_{t\mid t-1}$ in the EKF is time-varying. In the remainder, we will denote this constant $(1-\gamma)\beta$ by $\alpha$ for notational brevity.

Our complementary filter implementation is summarized in Algorithm 5. Note again the similarity to Algorithm 3. To stress this similarity even more, we call the update (4.60) the time update of the complementary filter and (4.4) the measurement update.

4.5 Evaluation based on experimental and simulated data

In this section, we apply the algorithms described in §4.1–§4.4 to both simulated and experimental data. Some general characteristics of the orientation estimation algorithms will be illustrated and the quality of the different algorithms will be analyzed. The simulated data allows for controlled analysis of the workings of the algorithms. Furthermore, it allows us to compare the different algorithms using Monte Carlo simulations. The experimental data shows the applicability to real-world scenarios. We will start by introducing the data sets.

Experimental data is collected using the setup shown in Figure 4.2, where data is collected using multiple mobile IMUs and smartphones. The algorithms presented in this section can be applied to measurements from any of these devices. However, we focus our analysis on the data from the Trivisio Colibri Wireless IMU [148]. In Figure 4.3, the inertial and magnetometer measurements from this IMU are displayed for around $100$ seconds during which the sensor is rotated around all three axes. The experiments are performed in a lab equipped with multiple cameras [151], able to track the optical markers shown in Figure 4.2. This provides highly accurate ground truth reference position and orientation information, against which we can compare our estimates. For comparison, the optical and IMU data need to be time-synchronized and aligned. We synchronize the data by correlating the norms of the gyroscope measurements and of the angular velocity estimated by the optical system. Alignment is done using the orientation estimates in combination with Theorem 4.2 from [59].

In Figure 4.4, simulated inertial and magnetometer measurements are displayed. The data represents a sensor that is kept stationary for 100 samples, after which it is rotated around all three axes. The sensor is assumed to be rotated around the origin of the accelerometer triad. Hence, during the entire data set, the accelerometer is assumed to only measure the gravity vector. The magnitude of the simulated gravity vector is $9.82~{}\metrepersquaresecond$. The magnitude of the simulated local magnetic field is equal to one. Its direction is approximately equal to that in Linköping, Sweden, where a dip angle of $71^{\circ}$ leads to a magnetic field $m^{\text{n}}=\begin{pmatrix}0.33&0&-0.95\end{pmatrix}^{\mathsf{T}}$. The simulated noise levels are

[TABLE]

Note that we deliberately chose the noise levels to be fairly high, to clearly illustrate the workings of the different algorithms.

Although our algorithms parametrize the orientation as quaternions, it is typically more intuitive to visualize the orientation estimates in Euler angles. Hence, we visualize our results in terms of roll, pitch and heading (yaw) angles. Both for the experimental data and for the simulated data, we are able to compare our estimates $\hat{q}^{\text{nb}}_{t}$ to reference orientations denoted $q^{\text{nb}}_{\text{ref},t}$. To represent the orientation error, we compute a difference quaternion $\Delta q_{t}$ as

[TABLE]

which can be converted to Euler angles for visualization. Note that using this definition, the orientation errors in Euler angles can be interpreted as the errors in roll, pitch and heading.

4.5.1 General characteristics

In this section, we will discuss some general characteristics of the orientation estimation problem and illustrate them in three different examples. Our goal is not to compare the different estimation algorithms, but to illustrate some characteristics common to all of them.

In Example 4.1 we focus on the accuracy of the orientation estimates that can be obtained if the state space model (3.75) is completely true. We illustrate that it is typically easier to obtain accurate roll and pitch estimates than it is to obtain accurate heading estimates.

Example 4.1** (Orientation estimation using inertial and magnetometer data)**

The orientation errors from the smoothing optimization approach in Algorithm 1 using simulated inertial and magnetometer measurements as illustrated in Figure 4.4 are depicted in the top plot of Figure 4.5. For comparison we also show the orientation errors from dead-reckoning the gyroscope measurements in the bottom plot (see also §1.2). These errors can be seen to drift over time.

Although the accelerometer and the magnetometer measurement noises are of equal magnitude, the heading angle is estimated with less accuracy compared to the roll and pitch angles. The reason for this is twofold. First, the signal to noise ratio for the magnetometer is worse than that of the accelerometer, since the magnetometer signal has a magnitude of $1$ while the accelerometer signal has a magnitude of $9.82~{}\metrepersquaresecond$. Second, only the horizontal component of the local magnetic field vector provides heading information. This component is fairly small due to the large dip angle ($71^{\circ}$) in Linköping, Sweden.

**

The accelerometer provides inclination information, while the magnetometer provides heading information, see §3.4. In case only inertial measurements and no magnetometer measurements are available, the heading can only be estimated using the gyroscope measurements and will therefore drift over time. This is illustrated in Example 4.2.

Example 4.2** (Orientation estimation using only inertial measurements)**

The orientation errors from the smoothing optimization approach in Algorithm 1 using simulated inertial measurements as presented in Figure 4.4 can be found in Figure 4.6. The roll and pitch angles can be seen to be accurate, while the heading angle drifts significantly. To obtain the results in Figure 4.6, we used data with the same noise realization as in Example 4.1. Hence, we refer to Figure 4.5 for comparison to the orientation errors from dead-reckoning the gyroscope measurements. The drift in the heading angle can be seen to be similar to the drift from dead-reckoning of the gyroscope measurements.

**

The two examples above assume that our state space model (3.75) is an accurate description of the measurements. In practice, however, this is not always the case, for instance due to the presence of magnetic material in the vicinity of the sensor. In Example 4.3 we illustrate that if the state space model does not accurately describe the data, it is not possible to obtain accurate orientation estimates.

Example 4.3** (Orientation estimation in the presence of magnetic material)**

We simulate 400 samples of stationary data. Between samples 150 and 250, we simulate the presence of a magnetic material, causing a change in the magnetic field of $\begin{pmatrix}0.1&0.3&0.5\end{pmatrix}^{\mathsf{T}}$. In Figure 4.7, we show the adapted magnetometer data and the orientation estimates using the smoothing optimization approach from §4.1. As can be seen, the heading estimates show significant errors when the magnetic material is present. Depending on the amount of disturbance and the uncertainty in the inertial and magnetometer measurements, magnetic material can also result in errors in the roll and pitch estimates.

**

4.5.2 Representing uncertainty

So far, we have discussed the quality of the orientation estimates in three different examples. However, we did not discuss the uncertainty of the estimates. Algorithms 1–4 provide estimates of the uncertainties. We will now discuss how these uncertainties can be displayed and interpreted and highlight some difficulties with this.

Both the optimization and the EKF approaches discussed in §4.1–§4.3 compute the uncertainty of the estimates in terms of a covariance matrix. Let us use the more general notation $\operatorname{cov}(\hat{\eta}_{t}^{\text{n}})$ for the covariance of the orientation deviation states $\hat{\eta}_{t}^{\text{n}}$, $t=1,\ldots,N$ computed in Algorithms 1, 2 and 4, and $\operatorname{cov}(\hat{q}_{t}^{\text{nb}})$ for the covariance of the quaternion states $\hat{q}_{t}^{\text{nb}}$, $t=1,\ldots,N$ computed in Algorithm 3. If the states would be in normal, Euclidean space, the square root of the diagonal of these matrices would represent the standard deviation $\sigma$ of the estimates in the different directions. These could then be visualized by for instance plotting $3\sigma$ confidence bounds around the estimates.

One could imagine that an equivalent way of visualizing the orientation deviation uncertainties, would be to compute the $3\sigma$ bounds in terms of orientation deviations in each of the three directions as

[TABLE]

after which the bounds can be parametrized in terms of quaternions as

[TABLE]

The resulting estimates and bounds are visualized in terms of Euler angles in Figure 4.8 for simulated data similar to the data presented in Figure 4.4. The orientation and its covariance are estimated using Algorithm 4. As can be seen, the bounds are difficult to interpret due to the wrapping of the Euler angles.

As argued in [38], it is more intuitive to directly represent the uncertainty in terms of orientation deviations. The covariance $\operatorname{cov}(\hat{\eta}_{t}^{\text{n}})$ can be interpreted as the uncertainty in the roll, pitch and heading angles as illustrated in Example 4.4.

Example 4.4** (Orientation estimation using only inertial measurements (continued))**

Since the accelerometer provides only inclination information, in the case of Example 4.2 where magnetometer measurements are unavailable, we expect only the roll and pitch angles to be estimated with small uncertainty. In fact, we expect the uncertainty of the heading at $t=1$ to be equal to the uncertainty of the initial $\Sigma_{\text{i}}$ from §3.6 and to steadily grow over time, depending on the amount of gyroscope noise. In Figure 4.9, we plot the standard deviation $\sigma$ of the orientation estimates computed using the smoothing algorithm from §4.1 as the square root of the diagonal elements of $\operatorname{cov}(\hat{\eta}_{t}^{\text{n}})$. As can be seen, the standard deviation of the yaw angle at $t=1$ is indeed $20^{\circ}$ as modeled in §3.6. The increase in the uncertainty in the yaw angle exactly matches the increase of the uncertainty due to dead-reckoning.

**

From Example 4.4 it can be concluded that $\operatorname{cov}(\hat{\eta}_{t}^{\text{n}})$ seems to be an intuitive measure of the uncertainty of the orientation estimates. The covariance $\operatorname{cov}(\hat{q}_{t}^{\text{nb}})$ computed by Algorithm 3 relates to $\operatorname{cov}(\hat{\eta}_{t}^{\text{n}})$ as

[TABLE]

The relation (4.71) allows us to compare the orientation estimates and covariances from the different algorithms in more detail in Example 4.5.

Example 4.5** (Orientation estimation using only inertial measurements (continued))**

As discussed in Example 4.4, using only inertial measurements and no magnetometer measurements, we can only expect to be able to accurately estimate the inclination. The uncertainty of the heading estimates grows over time. We will now analyze the behavior of the different algorithms in more detail for this specific example. In Table 4.1, we show the root mean square error (RMSE) values over $100$ Monte Carlo simulations for Algorithms 1–4. In Figure 4.10, we also represent the orientation estimates from the four algorithms for one of these realizations. As can be seen from both Table 4.1 and Figure 4.10, as expected, the smoothing algorithm outperforms the other algorithms. However, more surprisingly, the EKF with quaternion states has much larger errors in the heading angle. In Figure 4.11, we also show the standard deviations of the estimates from all algorithms. As can be seen, the EKF with quaternion states over-estimates its confidence in the estimates of the heading direction. This can most likely be attributed to linearization issues.

**

From Example 4.5, it can be concluded that not properly estimating the covariance can have negative effects on the quality of the estimates. It can also be concluded from this section that covariances can best be represented in terms of $\operatorname{cov}(\hat{\eta}_{t}^{\text{n}})$ but that they are difficult to visualize in Euler angles. Because of that, in the remainder of this section, we will typically plot the uncertainty in a separate plot like in Figure 4.11.

4.5.3 Comparing the different algorithms

We use the experimental data presented in Figure 4.3 to assess the quality of the estimates from the different algorithms. We also have a short sequence of stationary data available that allows us to determine the gyroscope bias and the sensor noise covariances. The gyroscope bias is subtracted from the data to allow for better orientation estimates. The standard deviation of the gyroscope noise is experimentally determined to be $\sigma_{\omega}=4.9\cdot 10^{-3}$. As discussed in §3.4, the covariance matrices in the measurement models reflect not only the sensor noise, but also the model uncertainty. Both the assumptions that the sensor’s acceleration is zero and that the magnetic field is homogeneous are only approximately true. Because of this, we choose $\sigma_{\text{a}}=2.6\cdot 10^{-1}$ and $\sigma_{\text{m}}=2.5\cdot 10^{-1}$. These values are a factor 10 respectively 100 larger than the experimentally determined noise values. Note that this choice is highly dependent on the data set. The RMSE values as compared to the optical reference system for the different methods described in this chapter are summarized in Table 4.2. A good value for $\alpha$ in the complementary filter is experimentally determined to be $0.001$ for this specific data set. As an illustration of the estimates, the orientation estimates as obtained using the smoothing algorithm and the orientations from the optical reference system are shown in Figure 4.12.

It is of course difficult to draw quantitative conclusions based on only one data set. The RMSE values in Table 4.2 should therefore mainly be seen as an indication of the performance. Comparing the different algorithms amongst each other is hard. In fact, the algorithms perform differently with respect to each other for different choices of the covariance matrices. Because of this, we will study the accuracy of the different methods using simulated data instead.

We run 100 Monte Carlo simulations where the simulated data illustrated in Figure 4.4 is generated with different noise realizations. Table 4.3 shows the mean RMSE for the five estimation algorithms. Algorithms 1–4 use the noise covariance matrices that are also used to generate the data. The complementary filter also uses these to determine the orientation from the accelerometer and magnetometer data. To combine this with the orientation from the gyroscope measurements, we need to determine a good value for $\alpha$. We empirically choose a value that leads to good performance for this data set. The orientation from the accelerometer and magnetometer measurements can be expected to be more accurate in the roll and the pitch than in the heading, see Example 4.1. Since $\alpha$ is scalar, however, it is not possible to weigh these contributions differently. Because of this, we present results both for $\alpha=0.07$, which leads to small RMSE values in the heading but relatively larger RMSE in the roll and the pitch, and for $\alpha=0.7$, which leads to small RMSE values in the roll and pitch but large RMSE in the heading.

From the results summarized in Table 4.3, it can be seen that the smoothing approach outperforms the filtering approaches. The estimates for one of the noise realizations are shown in Figure 4.13. For Algorithms 1–4, we also depict the covariance estimates in Figure 4.14. The filtering approaches from Algorithms 2–4 estimate the standard deviation of the orientation errors at $t=1$ to be equal to $20^{\circ}$. After this, they can be seen to converge to around $3.16^{\circ}$ degrees for the heading angle and $0.46^{\circ}$ for roll and pitch angles. The smoothing algorithm estimates an uncertainty in the heading angle of around $3.17^{\circ}$ for the first and last sample, while converging to a standard deviation of $2.25^{\circ}$ for the middle of the data set. For the roll and pitch angles, the initial and final uncertainties are estimated to be around $0.73^{\circ}$, converging to $0.39^{\circ}$ for the middle of the data set. Note that these values correspond fairly well with the RMSE values in Table 4.3.

For the Monte Carlo simulations described above, the three filtering algorithms perform similarly. However, differences can be seen when an update of the filter needs to correct the orientation estimates significantly. Examples for when this happens are when the initial orientation is not accurately known or when magnetometer measurements are not available for a longer period of time. In these cases, the uncertainty of the state is large and large corrections to the state estimates are needed when measurements become available. To analyze this case in more detail, we assume that the estimate of the initial orientation $\breve{q}_{1}^{\text{nb}}$ is normal distributed around the true initial orientation with a standard deviation of $20^{\circ}$. Hence, we do not use the first accelerometer and magnetometer data for initialization. Note that a standard deviation of $20^{\circ}$ is equal to the uncertainty on the initial state assumed by the algorithms. The results for $100$ Monte Carlo simulations are summarized in Table 4.4. As can be seen, specifically the EKF with quaternion states and the complementary filter with $\alpha=0.07$ perform worse than Algorithm 2 and Algorithm 4 for this data.

Which algorithm to use is highly application-specific. However, in general it can be concluded that all five algorithms actually produce fairly good orientation estimates, assuming that the models from §3.7 are indeed valid. The smoothing algorithm performs better than the filtering approaches but it is also the most computationally expensive. The EKF with quaternion states and the complementary filter suffer from linearization issues when large orientation corrections need to be made or when magnetometer data is unavailable.

4.6 Extending to pose estimation

In §4.5, we have evaluated the performance of Algorithms 1–5 for orientation estimation. The estimation methods presented in §4.1–§4.3 can also be used to estimate the sensor’s pose using the state space model (3.74). Complementary filtering is predominantly used for orientation estimation and will therefore not be considered in this section.

The pose estimation problem can be written as a smoothing optimization problem as

[TABLE]

with $x_{t}=\begin{pmatrix}p_{t}^{\mathsf{T}}&v_{t}^{\mathsf{T}}&(\eta_{t}^{\text{n}})^{\mathsf{T}}\end{pmatrix}^{\mathsf{T}}$ and

[TABLE]

In this section, we will discuss some details about the workings of the pose estimation algorithm using this model. We will not go through a complete derivation of the four algorithms. However, the adaptations that are needed to use Algorithms 1–4 for pose estimation can be found in Appendix B.

An important observation is that $e_{\text{a,p},t}$ and $e_{\text{a,v},t}$ in (4.73d) and (4.73e) depend on the orientation $R_{t}^{\text{nb}}$. Because of this, the position, velocity and orientation states are coupled. The position measurements therefore do not only provide information about the position and velocity, but also about the orientation of the sensor. This is the reason why it is no longer essential to include magnetometer data and to assume that the acceleration is approximately zero. However, the accuracy of the orientation estimates depends on the movements of the sensor. This will be illustrated below. For this, we simulate 400 samples of inertial and position measurements for a non-rotating sensor with noise levels

[TABLE]

We consider four different sensor motions. The results in this section are based on the solution to the smoothing optimization problem (4.6). First, we simulate data assuming that the sensor is stationary. For this case, the position measurements provide information about the inclination of the sensor, but not about its heading. This is illustrated in Example 4.6.

Example 4.6** (Pose estimation for a stationary sensor)**

We estimate the pose of a stationary sensor using simulated data and a smoothing algorithm that solves (4.6) as described in §4.1. The orientation error for a specific noise realization is depicted in Figure 4.15(a). The inclination errors can be seen to be small, while the heading estimates drift.

**

Next, in Example 4.7, we consider the case where the sensor has a constant linear acceleration. For this case, a drift in the orientation estimates can be seen in the direction that is orthogonal to the direction of the accelerometer measurements.

Example 4.7** (Pose estimation for a sensor with constant linear acceleration)**

We estimate the pose of a sensor with an acceleration of $1~{}\metrepersquaresecond$ in the $y$-direction using simulated data and obtain smoothed estimates by solving (4.6). The orientation error for a specific noise realization is depicted in Figure 4.15(b). Again, a drift can be seen in the orientation estimates. This drift is no longer only in the heading direction, but there is also a small drift on the roll angle.

**

Finally, in Example 4.8 we consider the case of a time-varying linear acceleration. Based on simulated data, we show that accurate heading estimates can be obtained for this case. Furthermore, we show that the larger the acceleration, the more accurate the heading estimates will be.

Example 4.8** (Pose estimation for a sensor with time-varying linear acceleration)**

We estimate the pose of a sensor with an acceleration in the $y$-direction of $\left(y_{\text{a},t}\right)_{y}\sim\mathcal{N}(0,0.5)~{}\metrepersquaresecond$ using simulated data and compute smoothed estimates by solving (4.6). The orientation error for a specific noise realization is depicted in Figure 4.15(c). Furthermore, we simulate data with $\left(y_{\text{a},t}\right)_{y}\sim\mathcal{N}(0,5)~{}\metrepersquaresecond$. The orientation errors based on this data can be found in Figure 4.15(d). As can be seen, for these cases, it is possible obtain reliable heading estimates using the state space model (3.74). The larger the acceleration, the more accurate the heading estimates.

**

In general, it can be concluded that it is possible to estimate position and orientation using the state space model (3.74). Except in the cases of constant or zero acceleration, it is possible to obtain drift-free orientation estimates. The heading accuracy depends on the amount of acceleration. This is summarized in Table 4.5 where the mean RMSE of the state estimates over 100 Monte Carlo simulations is shown. Four cases were considered, inspired by Examples 4.6–4.8.

Chapter 5 Calibration

In Chapter 4, we assumed that the sensors were properly calibrated. In practice, however, there are often calibration parameters to be taken into account. Examples of calibration parameters are the inertial sensor biases discussed in Chapter 2. Furthermore, calibration is specifically of concern when combining the inertial data with other sensors. In these cases, it is important that the inertial sensor axes and the axes of the additional sensors are aligned. Examples include using inertial sensors in combination with magnetometers [75, 126, 15] and with cameras [60, 83, 98].

In this chapter we will introduce several useful calibration methods. In §5.1 we explain how calibration parameters can be included as unknowns in the smoothing and filtering algorithms from §4.1–§4.3. This results in MAP estimates of the parameters. In §5.2 we instead focus on obtaining maximum likelihood (ML) estimates of the parameters. In §5.3, the workings of the calibration algorithms are illustrated by considering the gyroscope bias to be unknown in the orientation estimation problem. Finally, in §5.4, we discuss the topic of identifiability. Parameters are said to be identifiable if they can be estimated from the available data.

5.1 Maximum a posteriori calibration

As discussed in §3.1, unknown parameters $\theta$ can be estimated in the smoothing problem (3.1) as

[TABLE]

with

[TABLE]

Recall that a discussion on the choice of the prior of the parameters $p(\theta)$ and the states $p(x_{1})$ can be found in §3.6.

Within the filtering context we typically model the parameters as slowly time-varying states. These can be estimated by solving

[TABLE]

where

[TABLE]

Note that compared to §3.1, in (5.4) we do not consider the parameters to be part of $x_{t}$ but instead represent them explicitly. A prior $p(\theta_{1})$ on the parameters at $t=1$ has to be included as well as a dynamic model of the parameters.

Both the formulations (5.1) and (5.3) compute MAP estimates of the parameters. Algorithms 1–4 presented in Chapter 4 can straightforwardly be extended to also estimate these unknown parameters $\theta$ or $\theta_{1:N}$. This is illustrated in Example 5.1 for the case of orientation estimation in the presence of an unknown gyroscope bias. Note that Algorithm 5 can also be extended to estimate an unknown gyroscope bias. We will, however, not consider this. Instead, we refer the reader to [90].

Example 5.1** (MAP estimates of the gyroscope bias)**

It is possible to estimate an unknown gyroscope bias in the state space model (3.75). For this, the dynamic model (3.75a) in the smoothing problem described in §4.1 is assumed to include a constant gyroscope bias $\delta_{\omega}$ as

[TABLE]

The smoothing algorithm presented in §4.1 can be extended to also estimate $\delta_{\omega}$. Furthermore, the filtering algorithms presented in §4.2 and §4.3 can be extended to estimate $\delta_{\omega,t}$ for $t=1,\ldots,N$. Only minor changes to the algorithms presented in these sections are needed. These mainly concern including derivatives with respect to the additional unknowns $\delta_{\omega}$ or $\delta_{\omega,t}$. Explicit expressions for these can be found in Appendix C.

**

5.2 Maximum likelihood calibration

Alternatively, it is possible to obtain ML estimates of the parameters $\theta$ as

[TABLE]

Here, $\Theta\subseteq\mathbb{R}^{n_{\theta}}$ and $\mathcal{L}(\theta;y_{1:N})$ is referred to as the likelihood function. It is defined as $\mathcal{L}(\theta;y_{1:N})\triangleq p_{\theta}(Y_{1:N}=y_{1:N})$, where $Y_{1:N}$ are random variables and $y_{1:N}$ are a particular realization of $Y_{1:N}$. Using conditional probabilities and the fact that the logarithm is a monotonic function we have the following equivalent formulation of (5.6),

[TABLE]

where we use the convention that $y_{1:0}\triangleq\emptyset$. The ML estimator (5.7) enjoys well-understood theoretical properties including strong consistency, asymptotic normality, and asymptotic efficiency [82].

Due to the nonlinear nature of the orientation parametrization, our estimation problems are nonlinear, implying that there is no closed form solution available for the one step ahead predictor $p_{\theta}(Y_{t}=y_{t}\mid Y_{1:t-1}=y_{1:t-1})$ in (5.7). However, similar to the filtering approaches from Chapter 4, it is possible to approximate the one step ahead predictor according to

[TABLE]

where $\hat{y}_{t\mid t-1}(\theta)$ and $S_{t}(\theta)$ are defined in (4.39) and (4.38), respectively. Inserting (5.8) into (5.7) and neglecting all constants not depending on $\theta$ results in the following optimization problem,

[TABLE]

Unlike the optimization problems discussed so far, it is not straightforward to obtain an analytical expression of the gradient of (5.9). This is because it is defined recursively through the filtering update equations. In [3] and [134], different approaches to derive analytical expressions for objective functions of the same type as (5.9) are provided. They, however, consider the case of a linear model. Some methods for obtaining ML estimates of parameters in nonlinear models are explained in the tutorial by [131].

Instead of deriving analytical expressions for the gradient of (5.9), it is also possible to compute a numerical approximation of the gradient. Numerical gradients can be used in a number of different optimization algorithms, such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method, see e.g. [103]. Similar to the Gauss-Newton method, in the BFGS method, the parameters are iteratively updated until convergence. However, instead of using the Hessian approximation (4.8), BFGS iteratively estimates the Hessian using information from previous iterations. Hence, solving (5.9) using BFGS with numerical gradients, requires running at least $n_{\theta}+1$ filtering algorithms for each iteration. These are required to evaluate the objective function and to compute the numerical gradients. More evaluations can be necessary to compute a step length, see also §4.1.

Example 5.2** (ML estimates of the gyroscope bias)**

To obtain ML estimates of the gyroscope bias, we run the EKF with orientation deviation states from Algorithm 4 to obtain $\hat{y}_{t\mid t-1}(\delta_{\omega})$ and $S_{t}(\delta_{\omega})$ for a given value of $\delta_{\omega}$. This allows us to evaluate the objective function in (5.9). To compute $\hat{\delta}_{\omega}$, the optimization problem (5.9) is solved iteratively using BFGS.

**

5.3 Orientation estimation with an unknown gyroscope bias

We estimate the gyroscope bias in simulated data as described in §4.5 and illustrated in Figure 4.4. Compared to the data presented in §4.5, however, a constant gyroscope bias to be estimated is added. Using Monte Carlo simulations of this data, we illustrate a few specific features of the different ways to estimate the bias.

First, we focus on obtaining MAP estimates of the bias using the smoothing and filtering approaches as described in §5.1. We simulate the measurement noise as described in §4.5 and simulate the gyroscope bias as

[TABLE]

Note that $\sigma_{\delta_{\omega}}$ is a factor 10 larger than the value discussed in §3.6 to clearly illustrate the effect of the presence of a gyroscope bias. The priors $p(\theta)$ and $p(\theta_{1})$ in the smoothing and filtering algorithms are set equal to the distribution in (5.10). The covariance of the random walk model (5.5c) is set as $\Sigma_{\delta_{\omega,t}}=\sigma_{\delta_{\omega,t}}^{2}\mathcal{I}_{3}$ with $\sigma_{\delta_{\omega,t}}=1\cdot 10^{-10}$. This small value ensures that after convergence, the bias estimate is quite constant. The resulting mean RMSEs of the orientation over 100 Monte Carlo simulations are summarized in Table 5.1. Since the filtering algorithms typically have similar characteristics as discussed in §4.5, we only consider the EKF with orientation deviation states here. Comparing these results to the ones presented in Table 4.3, the RMSEs of the smoothing optimization algorithm are almost the same as when there was no gyroscope bias present. However, the filtering results are worse. This is because the bias needs some time to be properly estimated. This is illustrated in Figure 5.1 where the gyroscope bias estimates and their uncertainties are shown for the filtering algorithm.

A major difference between the MAP and the ML approaches, is that the MAP takes into account a prior on the gyroscope bias. We analyze the effect of this prior using $500$ Monte Carlo simulations, simulating the gyroscope bias to be $\delta_{\omega}=\begin{pmatrix}0.05&0.01&-0.04\end{pmatrix}^{\mathsf{T}}~{}\radianpersecond$. We study the estimated gyroscope biases using ML estimation, and using MAP estimation by including the gyroscope bias as an unknown in the smoothing algorithm. The smoothing algorithm assumes two different priors on the gyroscope bias $\delta_{\omega}\sim\mathcal{N}(0,\sigma^{2}_{\delta_{\omega}}\,\mathcal{I}_{3})$. In the first case, the prior on the gyroscope bias can well describe the data ($\sigma_{\delta_{\omega}}=0.05$). In the other case, the prior is too tight ($\sigma_{\delta_{\omega}}=1\cdot 10^{-3}$). The mean and standard deviations for the gyroscope bias estimates are summarized in Table 5.3. As can be seen, when the prior is chosen appropriately, the ML and MAP estimates are comparable. If the prior is too tight, the MAP estimates can be seen to be biased towards zero.