Kalman Filtering for Attitude Estimation with Quaternions and Concepts from Manifold Theory.

The problem of attitude estimation is broadly addressed using the Kalman filter formalism and unit quaternions to represent attitudes. This paper is also included in this framework, but introduces a new viewpoint from which the notions of "multiplicative update" and "covariance correction step" are conceived in a natural way. Concepts from manifold theory are used to define the moments of a distribution in a manifold. In particular, the mean and the covariance matrix of a distribution of unit quaternions are defined. Non-linear versions of the Kalman filter are developed applying these definitions. A simulation is designed to test the accuracy of the developed algorithms. The results of the simulation are analyzed and the best attitude estimator is selected according to the adopted performance metric.

1. Introduction

Mechanical state estimation of a vehicle is a field of interest. A vehicle is considered a rigid body, and its state of motion is represented by 4 mathematical objects: two of them represent its position and velocity, and the other two represent its orientation, and angular velocity. This paper is focused on the estimation of the angular state, composed of orientation, and angular velocity.

Although there are other mathematical tools used for estimation [1], the Kalman Filter [2] has become the algorithm par excellence in this area. Because of its simplicity, the rigor and elegance in its mathematical derivation, and its recursive nature it is very attractive for many practical applications. Its non-linear versions have been widely used in orientation estimation: the Extended Kalman Filter (EKF), and the Unscented Kalman Filter (UKF) [3]. However, there are problems arising from the used parametrization to represent the orientation.

The orientation of a system is represented by the rotation transformation that relates two reference frames: a reference frame anchored to that system, and an external reference frame. A thorough survey of attitude representations is provided in Reference [4]. The parametrization used to represent the rotation transformation could be singular, or present discontinuities among others. Table 1 summarizes the main characteristics of the most used parametrizations.

Table 1

Main characteristics of the most used parametrizations to represent an orientation.

Representation	Parameters	Continuous	Non-Singular	Linear Evolution Equation
Euler angles	3	✗	✗	✗
Axis-angle	3–4	✗	✗	✗
Rotation matrix	9	✓	✓	✓
Unit quaternion	4	✓	✓	✓

Having in mind that the special orthogonal group SO(3) has dimension three, ideally we would seek for a continuous and non-singular representation expressed by 3 parameters. However, since 1964 we know that “...it is topologically impossible to have a global 3-dimensional parametrization without singular points for the rotation group” [5]. Knowing this, we would not be wrong to say that unit quaternions are the most convenient representation we have, and that we will have for orientations. In Reference [6] the literature on attitude estimation is reviewed until 1982, when other parametrizations like Euler angles were common, and founds the basis of modern quaternion-based attitude estimation, in which this paper is supported. After that work, many others have explored this viewpoint, and have demonstrated its superiority [7,8,9,10,11,12].

Quaternions are 4-dimensional entities, but only those having unit norm represent a rotation transformation. This fact implies a problem in applying the ordinary Kalman Filter, so different approaches have emerged. Since a quaternion is of dimension 4, one tends to think at first on a covariance matrix, and in the direct application of the Kalman Filter [13]. Given that all predictions are contained in the surface defined by the unit constraint, the covariance matrix shrinks in the orthogonal direction to this surface, which leads to a singular covariance matrix after several updates. A second perspective was firstly approached in Reference [6] and was after named as “Multiplicative Extended Kalman Filter” [8,11,12]. In this second approach we define an “error-quaternion” that is transformed to a 3-vector. We use this vector to build the covariance matrix, and we talk about a “ representation of the quaternion covariance matrix”. However, there are still details in this adaptation that are currently being developed. Namely, the “covariance correction step” [14].

This paper presents a new viewpoint on the problem of attitude estimation using Kalman filters when the orientation is represented by unit quaternions. Noticing that unit quaternions live in a manifold (the unit sphere in ), we use basic concepts from manifold theory to define the mean and covariance matrix of a distribution of unit quaternions. With these definitions we develop two estimators based on the Kalman filter (one EKF-based and another UKF-based) arriving at the concepts of “multiplicative update” and “covariance correction step” in a natural and satisfying way. The inartificial emergence of these ideas establishes a solid foundation for the development of general navigation algorithms. Lastly, we also analyze the accuracy in the estimations of these two estimators using simulations.

The organization of this paper is as follows. In Section 2 we review quaternion basics. We also expose the new viewpoint on the definition of the quaternion mean and covariance matrix. In Section 3 we present the developed estimation algorithms. In Section 4 we define the performance metric, describe the simulation scheme, and present the results of the simulations. We also discuss the results. Finally, Section 5 concludes the paper.

2. Quaternions Describing Orientations

2.1. Quaternions

Quaternions are hypercomplex numbers composed of a real part and an imaginary part. The imaginary part is expressed using three different imaginary units satisfying the Hamilton axiom:

A quaternion can be represented with 4 real numbers, and using several notations:

We will denote quaternions with bold italic symbols (), while vectors will be denoted with bold upright symbols (). Vectors will be written in matrix form, and the transposed of a matrix will be denoted as .

Quaternion product is defined by Equation (1) which produces the multiplication rule where represents the usual dot product, and represents the 3-vector cross product. Note that the quaternion product (*) is different from the product denoted by (⊗) in reference [4]. Given this multiplication rule, the inverse of a quaternion (the one for which ) is given by where represents the complex conjugate quaternion. Note that if is a unit quaternion (a quaternion with ), then .

2.2. Quaternions Representing Rotations

Each rotation transformation is mapped with a rotation matrix and with two unit quaternions and all of them related through

Note that .

Quaternions representing rotations have the form where denotes the unit vector that defines the rotation axis, and the angle of rotation. Having this form, they satisfy the restriction

This means that quaternions describing rotations live in the unit sphere of , . This space is a manifold, and some concepts regarding these mathematical objects are useful in our context. In particular, the concept of chart is of special interest.

2.3. Distributions of Unit Quaternions

When dealing with the Kalman filter, the distribution of a random variable is encoded by its mean and its covariance matrix defined as

This definition makes sense when our random variables are defined in the Euclidean space. But how do we define the covariance matrix of a random variable living in a manifold like ours? How can we define the covariance for unit quaternions if does not represent a rotation? (Unit quaternions form a group under multiplication, but not under addition. This means that the addition of two unit quaternions may not result in another unit quaternion. Therefore, the addition of two unit quaternions may not represent a rotation.) What would be the covariance matrix if each quaternion was equiprobable in the unit sphere? We cannot redefine the covariance matrix, because the Kalman filter uses this precise form in its derivations, but we can take advantage of the properties of a manifold. Let us retrieve some important definitions:

A homeomorphism is a function If such a function exists, we say that X and Y are  homeomorphic.

f is a bijection,

f is continuous,

its inverse function is continuous.

A n-manifold

A chart for a manifold with φ a homeomorphism. Traditionally, a chart is expressed as the pair

Given these definitions we can continue our reasoning.

In Reference [8] it talks about four “attitude error representations”. Namely, the one we will call Orthographic (O), the Rodrigues Parameters (RP), the Modified Rodrigues Parameters (MRP), and the Rotation Vector (RV). The first three are what we know as stereographic projections (and are called Orthographic, Gnomonic, and Stereographic respectively). The last one is a projection called Equidistant. But all four are charts defining a homeomorphism from the manifold to the Euclidean space . This is, they map a point in the manifold with a point in . Table 2 arranges these chart definitions, together with their domain and image. We must ensure the charts to be bijections so that they properly define a homeomorphism, and that they do not map and with different points of since they represent the same rotation. We achieve this by the given definitions of the domain and image for each chart.

Table 2

Main characteristics of the charts studied.

Chart	Domain	Image	e=φ(q)	q=φ^−1(e)
O	{q∈S^3:q_0≥0}	{e∈R^3:∥e∥≤2}	2q	1−∥e∥^24e/2
RP	{q∈S^3:q_0>0}	R^3	2qq_0	14+∥e∥^22e
MRP	{q∈S^3:q_0≥0}	{e∈R^3:∥e∥≤4}	4q1+q_0	116+∥e∥^216−∥e∥^28e
RV	{q∈S^3:q_0≥0}	{e∈R^3:∥e∥≤π}	2q^arcsin∥q∥	cos∥e∥2e^sin∥e∥2

Figure 1 shows how points in the sphere (subspace of the sphere where quaternions live) are mapped to points in (subspace of where the images of the charts are contained) through each one of the named charts. Since our charts are homeomorphisms, it is possible to invert the functions. Figure 2 shows how points from are mapped to points in the manifold through the inverted charts. As pointed in Reference [8], all four charts provide the same second-order approximation for a point near the origin, to a quaternion

Figure 1

Points in the manifold with are mapped with points in the Euclidean space through each chart . (a) ; (b) O; (c) RP; (d) MRP; (e) RV.

Figure 2

Points in the Euclidean space with are mapped with points in the manifold through each chart inverse . (a) ; (b) O; (c) RP; (d) MRP; (e) RV.

We should notice that having and different metrics, a chart will inevitably produce a deformation of the space. However, for quaternions in the neighborhood of the identity quaternion (top of the sphere), our charts behave like the identity transformation between the imaginary part of these quaternions, and the points near the origin in as suggested by (10). This is a desirable property, as this means that the space around the identity quaternion closely resembles the Euclidean space, which is the space for which the Kalman filter is designed. But this just happens in the neighborhood of the identity quaternion. However, we can extend this property for any quaternion noting that any quaternion can be expressed as a “deviation” from the first one through the quaternion product: where represents such a deviation. (This definition is arbitrary: we could have chosen to relate the quaternions through , but it is important to establish one of these definitions, and then be consequent with it. However, (11) entails a computational advantage for the computation of (37).) Then, we define a chart for each quaternion as where and where we have denoted the point of the Euclidean space mapped with the quaternion through the chart as . Then, we will have a set of charts , each one resembling the Euclidean space around a quaternion , and mapping this last quaternion to the origin of . We will refer to the Euclidean space associated with the chart as the -centered chart. Thus, the homeomorphism takes a point in the -centered chart and maps it to a point in the manifold through

After reviewing these concepts, we can define the covariance matrix of a distribution of unit quaternions.

Given a unit quaternion and a chart , we will define the expected value of a distribution of unit quaternions in the -centered chart as and its covariance matrix as and the probability density of each unit quaternion would be defined through the homeomorphism . Then, a distribution of unit quaternions needs of four mathematical objects to be encoded: . Although a distribution of unit quaternions is unique, given this definition, its expected value and its covariance matrix may take different values depending on the chosen quaternion and chart . However, knowing that the Kalman filter is designed for the Euclidean space, it will be convenient to choose a unit quaternion central in the distribution, in order that the manifold space around it closely resembles the most significant region for the covariance matrix in the -centered chart. It is particularly convenient to choose a quaternion such that so that the covariance matrix is centered in the origin of the -centered chart.

2.4. Transition Maps

At some step of the Kalman filter, we will have a distribution of unit quaternions defined in a -centered chart, and we will be interested in expressing our distribution in another -centered chart. The concept of transition map is relevant for this purpose.

Given two charts with transition map. Being that

For the present case, let us consider two unit quaternions and both related through

These two quaternions define the charts and . We build the transition map that relates a point expressed in the -centered chart with a point expressed in the -centered chart doing

That is to say, first we take the point in the -centered chart, and we obtain its associated quaternion in the manifold using . Then, we transform this quaternion to a point in the -centered chart. Nevertheless, knowing the quaternion we do not need to explicitly compute . In fact, being able to express the same quaternion as two different deviations,

Note the equivalence of expressions (18c) and (19).

Table 3 displays the transition maps for the charts studied. The detailed derivations of these transition maps can be found in Appendix A. Figure 3 attempts to provide some insight into how points are transformed through the transition map of each chart.

Table 3

Transition maps for the charts studied.

Chart	Transition Map e^p¯e^q¯
O	δ¯_0e^q¯−4−∥e^q¯∥^2δ¯−δ¯×e^q¯
RP	2δ¯_0e^q¯−2δ¯−δ¯×e^q¯2δ¯_0+δ¯·e^q¯
MRP	48δ¯_0e^q¯−(16−∥e^q¯∥^2)δ¯−8δ¯×e^q¯16+∥e^q¯∥^2+δ¯_0(16−∥e^q¯∥^2)+8δ¯·e^q¯
RV	2δ^p¯∥δ^p¯∥arcsin∥δ^p¯∥,   with δ^p¯=δ¯_0e^^q¯sin∥e^q¯∥2−cos∥e^q¯∥2δ¯−δ¯×e^^q¯sin∥e^q¯∥2

Figure 3

Points in a -centered chart are transformed using the transition map defined by each chart, and travel to the chart centered in the quaternion mapped with in the previous chart. (a) ; (b) O; (c) RP; (d) MRP; (e) RV.

3. Manifold Kalman Filters

In this section we present the models adopted for the Manifold Kalman Filters (MKF), and we display the resulting algorithms.

The state of the system at a time t is defined by an orientation, encoded with a unit quaternion and by an angular velocity . We will consider them to be random variables, and we will try to estimate their value using a Kalman filter.

Our unit quaternions will define the rotation transformation that relates a vector expressed in a reference frame attached to the solid whose state we want to describe, with the same vector expressed in an external reference frame

For example, if we measure an acceleration in reference frame the acceleration in the inertial reference frame would be given by . This acceleration would be the one that we would have to integrate to obtain the position estimated by an accelerometer.

The vector will define the angular velocity of the solid measured in . Note that we do not include the bias of the sensors in the state of our system. We will assume that our sensors are calibrated, so the biases are zero.

We can predict the value of the random variables that describe the state of our system through the following motion equations: where is a random variable that represents the process noise, and is associated with the torque acting on the system, and with its inertia tensor. Its expected value at a given time t will be denoted as and its covariance matrix will be denoted as .

We will assume that we have sensors giving measurements of angular velocity (which provide information about the relative change in orientation), and of a vector whose value expressed in the external reference frame is known (this provides information about absolute orientation). Examples of such sensors could be a gyroscope giving angular velocity measurements, an accelerometer measuring the gravity vector near the Earth surface (), or a magnetometer measuring the Earth magnetic field (). The measurement model relates these measurements with the variables that describe the state of the system: where and are random variables with zero mean and covariance matrices and respectively that represent the measurement noises, and is another random variable with mean and covariance matrix representing external disturbances in the measurement of the vector . For example, it could represent accelerations others than gravity for an accelerometer, or magnetic disturbances produced by moving irons for a magnetometer.

We will assume that the measurements arrive at discrete times . The format will be used to denote a variable x at a time t, having included measurements up to a time with . For the n-th time stamp, in which a measurement arrives, we will write for the sake of simplicity. Then, our knowledge about the state at a time t, having included measurements up to a time with is described by a distribution encoded in the collection of mathematical objects as described in Section 2.3. For the present case, is the expected value of the distribution, and is its covariance matrix, both expressing the quaternion distribution in the -centered chart. Preferably, will be a unit quaternion central in the distribution, so that the mapping of points from the -centered chart to the manifold causes minimal deformation in such distribution. The unit quaternion will be our best estimation of the real quaternion that defines the orientation of the system with respect to the external reference frame at time t.

The following subsections present the developed Kalman filters: one version based on the EKF and another version based on the UKF. The EKF is based on the linearization of the non-linear models to calculate the predicted covariance matrices. That is, the EKF approximates non-linear functions using their Jacobian matrices. To apply the EKF, our functions must be differentiable. On the other hand, the UKF is based on a deterministic sampling to approximate the distribution of our random variables. We select a minimal set of samples whose mean and covariance matrix are those of the state distribution. Then, they are transformed by the non-linear models, and the resulting set of points is used to compute the means and covariance matrices necessary to perform the Kalman update. This second approach does not need the functions to be differentiable.

3.1. Manifold Extended Kalman Filter

In this section we present the EKF-based estimator: the Manifold Extended Kalman Filter (MEKF). We offer here the main results of the more detailed derivation given in Appendix B.

A measurement arrives at time . Our knowledge about the orientation at a previous time is described by a distribution expressed in the -centered chart. We assume that this distribution has mean and covariance matrix . This is, we have an initial four

The state prediction at time given all the information up to is computed through with

The measurement prediction at the same time is given by where stands for

At this point, we compute the Kalman gain and use it to obtain the optimal estimation of the state: where . Finally, we need to obtain the updated unit quaternion, and compute the mean and the covariance matrix in the -centered chart, so that the distribution is expressed in the same conditions as at the beginning of the iteration. The point that results from (41), and that is defined in the -centered chart, correspond to a unit quaternion in the manifold. This is the updated unit quaternion which we are looking for:

Knowing that the Kalman update (41) could produce any point in the -centered chart we will need to “saturate” to the closest point contained in the image of each chart. The point in the -centered chart is the origin in the -centered chart. Then, the expected value of the state in this new chart will be given by as at the beginning of the iteration.

To update the covariance matrix we need to consider its definition (15). We want to compute having and knowing the relation provided by the transition maps in Table 3. Continuing with the EKF philosophy, the update for the covariance matrix will be found by linearizing around the point where the majority of information is comprised (in our case, the point ): where we have used the big O notation to describe the limiting behavior of the error term of the approximation as . In particular, if we define then, and the final update for the covariance matrix will be computed through

Table 4 summarizes the resulting -matrix for each chart, along with their application domain. A detailed derivation of these -matrices can be found in Appendix C.

Table 4

-matrices for the transition maps of the charts studied.

Chart	T(δ¯) Matrix	Domain
O	δ¯_0I−δ¯_×+δ¯δ¯^Tδ¯_0	{δ¯∈S^3:δ¯_0>0}
RP	δ¯_0δ¯_0I−δ¯_×	{δ¯∈S^3:δ¯_0≠0}
MRP	121+δ¯_0δ¯_0I−δ¯_×+δ¯δ¯^T	{δ¯∈S^3:δ¯_0≥0}
RV	δ¯_0I−δ¯^δ¯^^T−δ¯_×∥δ¯∥arcsin∥δ¯∥+δ¯^δ¯^^T	{δ¯∈S^3:δ¯_0≥0,∥δ¯∥≠0}

After the final computation we obtain the four that is a condition equivalent to (27) in which we started the iteration.

3.2. Manifold Unscented Kalman Filter

In this section we present the UKF-based estimator: the Manifold Unscented Kalman Filter (MUKF).

A measurement arrives at time . Our knowledge about the orientation at a previous time is described by a distribution expressed in the -centered chart. This distribution is encoded in the four

The first step in the UKF is to create the augmented mean and covariance matrix . Since the measurement equations are linear for the random variables and we can leave their covariance matrices out of the augmented one and add them later:

Then, we obtain the matrix which satisfies and we use it to generate the sigma points as described in Ref. [15]: being for where regulates the importance given to the sigma point in the computation of the mean. These sigma points are expressed in the -centered chart. We need to express them in the manifold before applying the evolution equations and the measurement equations: where for the j-th sigma point, is its chart point part and is the quaternion with which it is mapped, is its angular velocity part, is its angular velocity noise part, is its angular velocity prediction, is the quaternion part of its prediction (we have assumed that the angular velocity is constant in the time interval so that we can use (A20)), is the vector process noise part, is its vector measurement prediction, is its angular velocity measurement prediction, and . Note that when applying the inverse chart we will need to “saturate” to the closest point in the image of . Having these new sigma points, we can obtain the means and covariance matrices of the distributions present in the UKF. First, defining the means are computed through where we have used a variation of the result provided in Ref. [16]. Namely, with for . This result is shown to minimize the fourth order approximation of the distance defined as the sum of squared angles between the rotation transformation represented by each quaternion and the one represented by . This approach to compute the mean quaternion is extremely efficient, and its derivation is elegant and simple. In order to ensure that it is useful to remember the property that both and represent the same rotation. This property is also useful for introducing the quaternions in the domain of to execute the next step of the filter.

After this, we use the obtained mean quaternion to express each sigma point in the -centered chart, and compute the covariance matrices: where we have denoted . Finally, we compute the UKF version of the Kalman gain and we use it to obtain the optimal estimation of the state: arriving at the same conditions in which we began the iteration, with a distribution expressed in the -centered chart, and encoded by the four

Our best estimation for the orientation at this time is being the part of the mean that represent the quaternion in the -centered chart.

Note that setting and at the beginning of each iteration yields the traditional version of the algorithm, where a “reset operation” is performed instead of the covariance matrix update.

4. Simulation Results

This section presents the results of the simulations used to measure the accuracy of each estimator. Simulations are chosen instead of real experiments because a real system entails an uncertainty in the measurement of the true attitude: the attitude that is used to compare with that estimated by the algorithms. There are sources of error ranging from a miscalibration of the measurement system to a possible bias in the “true attitude” produced by another attitude estimator, which makes it problematic to define an adequate metric to measure the accuracy of the algorithms. For this reason, the authors consider that using a simulation is more reliable to avoid possible biases in the results due to said sources of error. Others have performed similar types of tests [7,17]. However, the results do not seem to be statistically conclusive: only the estimations of some orientation trajectories are shown.

We perform our comparison through a simulation in which we do have an absolute knowledge of the attitude of the system: a true oracle exists in a simulation. Therefore, we can compare the real orientation with the attitude estimated by the algorithms having fed them only with simulated measurements that we obtain from such known orientations. We will extract our performance metrics from a wide set of orientation trajectories in order to obtain statistically conclusive results.

We try to answer three questions with the simulation test. The first question is, is there a chart for which we get a greater accuracy in attitude estimation? The second one is, what algorithm produces the most accurate attitude estimation, the MEKF or the MUKF? The last question stems from the fact that previous algorithms on attitude estimation, such as the Multiplicative Extended Kalman Filter, did not contemplate updating the distribution from one chart to another as done at (47b) in the MEKF. However, their estimators performed well [6,7,12]. Then the third question is, does this “chart update” imply an improvement in the accuracy of the attitude estimation?

Although a simulation has been used to compare our algorithms, these have also been tested with a real IMU. In the Supplementary Materials one can find a demonstration video, the source code used in the video, the source code used to generate the simulations, and the source code used to obtain the computational cost of the algorithms in each platform.

4.1. Performance Metric

We have already described a quaternion as a deviation from another quaternion as . Now we define the instantaneous error between an estimated attitude, represented by a unit quaternion and the real attitude, represented by the unit quaternion as the angle we have to rotate one of them to transform it into the other. This is, the angle of the rotation transformation defined by the quaternion such that . Recalling (6), this angle can be computed as: having previously ensured that using the fact that both and represent the same rotation transformation.

Angle will vary along an orientation trajectory. Then, we will define the mean error in orientation estimation for a given trajectory starting at time and ending at time as

Finally, will depend on the followed trajectory, and on the set of taken measurements. We will need to generate several orientation trajectories to obtain the mean value and the variance that characterize the distribution of the error in orientation estimation for each algorithm. We will define the confidence interval for the computed as where is the number of samples taken for the computation, so that is the variance of the sample mean distribution.

Being that the lower the better, the value of gives us a measure of how well an algorithm estimates the orientation. We will consider that the performance of an algorithm A is better than the performance of other algorithm B if and their confidence intervals do not overlap.

4.2. Simulation Scheme

To compute the performance metrics we will need to generate a large number of simulations. Each independent simulation will consist of three steps: initialization, convergence, and estimation.

In the initialization step we set up the initial conditions accordingly to the chosen simulation parameters. This includes generating the initial unit quaternion from a uniform distribution in setting the initial angular velocity to zero, setting the update frequency generating the variances of the process noises and from a uniform distribution in the intervals and respectively, and initializing the estimation algorithm. The initialization of the MEKF includes setting and . On the other hand, the initialization of the MUKF includes setting and . The angular velocity is not initialized to in the MUKF because it has been observed that it is sometimes necessary to “break the symmetry” for the algorithm to converge; especially when we do not apply the chart update (when we perform the “reset operation”) for the RV chart. The covariance matrices that appear in both algorithms are initialized as (“p.d.u.” stands for “Procedure Defined Unit”. In the present case it depends on the definition of the vector ), , , where and are the variances of the measurement noise that will be used in the simulation. We give this information about the measurement noise to the algorithms because it can be obtained offline, while the information about the process noise cannot. Given that a priori we cannot know how the system will behave, the values of and have been chosen according to what we understand could be normal. Choosing these values we are assuming that after a second it is normal for the angular velocity to have changed by and also that it is normal to find external noises added to the vector of magnitude . For the mean values we set and .

In the convergence step we keep the system in the initial orientation . Simulated measurements are generated using (23) and (24). For each measurement, a different is sampled from a uniform distribution in the unit sphere of . The values for each component of and are obtained from normal distributions with zero mean and variances and respectively. The term in (23) is obtained from the true attitude , which in the convergence step takes the value of . The term in (24) is the true angular velocity, which in the convergence step takes the value . The tested algorithm updates its state estimation until the inequality is satisfied, where is the value of the error (72), and is a parameter in the simulation. The convergence step could have been replaced by an initialization of the attitude estimated by the algorithm to the real value but then it would have also been necessary to fix a certain covariance matrix. Since the metric of the space generated by each chart is different, it is difficult to set a covariance matrix that provides the same information for each chart. It seemed more natural to the authors to allow the algorithm to find the true attitude by its own means, and for the covariance matrix to converge to a value in each case.

Finally, in the estimation step we generate a random but continuous orientation sequence using a Wiener process for the angular velocity: where is a random vector whose components are sampled from a normal distribution with zero mean and variance and is the simulation time step that is related to the algorithm time step Δt trough being an integer parameter that determines the simulation updates per algorithm update. Note that we multiply by and not by . We do it this way so that the covariance matrix after k steps does not depend on the simulation time step . In fact, after a time the covariance matrix of the angular velocity will have grown by and not by . After each simulation updates, a simulated measurement is generated in the same way it was done in the convergence step, and the algorithm is updated with it. The simulation will run for a time where is an integer number. This way we will perform the last algorithm update at the end of the simulation. The error (72) will be evaluated after each algorithm update, and it will be added up through the simulation to obtain the averaged error (73). After each simulation, we will obtain a sample for the computation of and . We will perform of these simulations to obtain the confidence interval (74).

4.3. Results

In this section we present the results of the simulations. The algorithms are tested for update frequencies in the interval . This range has been chosen thinking about the possible limitations of a real system. For example, the maximum data rate of a low cost IMU is around . On the other hand, the update frequency may be limited by processing. The computational cost of each estimator has been evaluated in two platforms: an Arduino MEGA 2560, and a Raspberry Pi 3 Model B. The code has been written in c++. The resulting maximum update frequencies are presented in Figure 4, which indicates that the MEKF can be executed approximately 3 times faster than the MUKF.

Figure 4

Maximum update frequency for each approach. The lines at the top represent the mean and the deviation () of the distribution of maximum update frequencies. “CU” stands for Chart Update, while “NCU” stands for No Chart Update.

Although the algorithms have been developed allowing a different for each update, the simulations are performed using a constant Δt, and the simulation parameters depicted in Table 5.

Table 5

Parameters used in the simulations.

Parameter	Value
θe0	1^∘
Tsim	10s
dtdtsim	100
N_s	1000
Qmaxω	10^2rads^2/s^3
Qmaxv	1p.d.u.
R	{10^−2,10^−4,10^−6}
R^ω	Rrads^2/s^2
R^v	Rp.d.u.
W_0	1/25

The parameters Tsim, dtdtsim, and have been chosen trying to reach a compromise between the precision of the results, and the execution time of the simulation. The values for and have been chosen in such a way that the estimation algorithms face both normal situations ( and ) and situations that were not foreseen ( or ). A typical low cost IMU has and . The values chosen for R represent an imprecise sensor (), a normal sensor (), and a precise sensor (). The value of has been chosen so that all sigma points have the same importance, but very similar results, if not identical, have been obtained for other selections of .

4.3.1. Chart Choice

The results of the simulation are presented in Figure 5. The average of the performance metric is shown along with its confidence interval for each of the selected update frequencies. The results of the MEKF and the MUKF are shown in different graphs, but drawn in the same one are the results for each chart and for a given MKF. In this way we are able to distinguish if a chart has an advantage over the others.

Figure 5

Mean of the performance metric for each approach. Results from different charts are plotted in the same graph. Results from different MKF are plotted in different graphs. Bars represent the confidence interval () for the mean computation.

We observe that there is no chart that is especially advantageous. All things being equal, we would opt for the RP chart. For this chart it is not necessary to worry about the domain since it maps and with the same point of and with the same -matrix; or of the image since it is all . In addition, the expressions of and the -matrix for the MEKF are simpler for the RP chart. These computational advantages make us prefer the RP chart over the others.

4.3.2. MEKF vs. MUKF

Figure 6 also presents the results of the simulations. This time, we display on the same graph the resulting performance metrics for the MUKF and the MEKF when the RP chart is used. In this way, we can distinguish if one MKF has an advantage over the other.

Figure 6

Mean of the performance metric for each MKF. Only results for the RP chart are plotted. Bars represent the confidence interval () in the mean computation.

We note that the MEKF performs the same or better than the MUKF. This differs from the usual experience, in which the UKF outperforms the EKF in traditional non-linear estimation applications. The fact that the charts resemble the Euclidean space near the origin (see Section 2.3) might be favoring the MEKF, since the Jacobian matrices, used to approximate the non-linear functions, are defined at that point. However, the sigma points generated for the MUKF are sampled far from the origin of the chart, where the non-linearities become notorious. We are facing a very particular scenario in which the model is approximately linear for the MEKF, while for the MUKF it is not. In addition, due to the difference in computational cost (see Figure 4), the MUKF update frequencies will generally be lower than those of the MEKF, which will imply worse accuracy in its estimations. Then, the MEKF with the RP chart seems to be our best option.

4.3.3. Chart Update vs. No Chart Update

Figure 7 presents the results of each MKF with each chart in a different graph, but displayed in the same one are the results using the “chart update” and the results without using it.

Figure 7

Mean of the performance metric for each approach. Results from the approach in which we apply the chart update, and those of which we do not apply it are plotted together. Bars represent the confidence interval () in the mean computation.

We can observe that there is almost no difference between using the “chart update” and not using it. The concepts used in this paper have helped us to understand the mechanisms of the MKF, and ultimately to arrive to the concepts of “multiplicative update”, and of “covariance correction step” with the -matrix definition. However, it is not necessary to apply the latest update (47b) in practice: we will obtain essentially the same accuracy in our estimations.

5. Conclusions

We have used concepts from manifold theory to define the expected value and the covariance matrix of a distribution in a manifold. In particular, we have defined the expected value and covariance matrix of a distribution of unit quaternions in , the unit sphere in , using the concept of chart. These definitions have helped us to develop Kalman filters for orientation estimation, where the attitude has been represented by a unit quaternion. They have also helped us solve the problem of the “covariance correction step”. Two estimators have been developed: one based on the EKF (the MEKF), and another based on the UKF (the MUKF). The MEKF and the MUKF have been tested in simulations, and some results have been obtained. The conclusions of the simulations are:

There is no chart that presents a clear advantage over the others, but the RP chart has some characteristics that motivate us to prefer it.

The MEKF is preferable to the MUKF due to its lower computational cost and its greater accuracy in orientation estimation.

The “chart update” is not necessary for the MKF in practice.

Then, the MEKF with the RP chart and without applying the “chart update” is our best attitude estimator according to the adopted performance metric. This algorithm resembles the conventional “Multiplicative Extended Kalman Filter”, but we have obtained the MEKF without having to redefine any aspect of the classic Kalman filter.