In this section, we describe the methods used for the proposed evaluation. First we detail the orientation estimation algorithm, which consists of a Kalman filter. Lastly, the procedure to emulate different magnetic field distributions from error-free turn rate measurements is described. In this subsection, we also detail how we add noise and biases to the error-free measurements in order to obtain similar quality measurements as if they were recorded with MEMS sensors.
2.1. Orientation Estimation Filter
The orientation estimation algorithm aims at fusing the information of gyroscopes, accelerometers and magnetometers in an optimal way to obtain the orientation of the sensor [
18]. The Kalman filter is suited for modelling continuous variables whose system model equations are linear and whose system and measurements noise are Gaussian [
19]. Taking into account the non-linearities of the equations involved in the orientation computation, it is convenient to use an extended Kalman filter [
6] or an unscented Kalman filter [
4,
17].
The state vector of the filter proposed in this work,
, is composed by the Euler angles and the biases of the gyroscopes
being
the roll angle,
the pitch angle,
the yaw angle, which are the Euler angles that represent the orientation, and
the biases of the gyroscopes, for
} representing the three orthogonal axes. Although the error state Kalman filter is more common in navigation applications, we have chosen the full state Kalman filter because both have similar performance, as demonstrated in [
20].
The diagram (
Figure 1) shows the computation taken for the time stamp
k. The variables
,
and
represent the measurements of the accelerometers, magnetometers and gyroscopes respectively. The states
as well as
are initialized to zero. The states
and
, as well as the initial value of
,
and
are computed during the calibration phase taking into account the average of 1
of static measurements. After the prediction stage, the state vector is represented as
. Next, some characteristics of the measurements are checked by the detectors, whose goal is to determine whether the updates can be applied or not. If the update stage takes place, the state vector is updated,
, and the
iteration finishes. If no update takes place, the predicted state vector
is outputted and the
iteration finishes.
2.1.1. Prediction Stage
The navigation frame we use in this work is the East-North-Up coordinate system (see
Figure 2). The gyroscopes measure in body frame the turn rate of the sensor with respect to the inertial frame,
. The outputted turn rate is composed of:
The term is the turn rate of the sensor measured in body frame with respect to the navigation frame. The term is the transport rate, which represents the turn rate of the navigation frame with respect to the Earth-fixed frame. The term represents the turn rate of the Earth with respect to the inertial frame. Lastly, the term is equivalent to , which represents the rotation of the navigation frame with respect to the inertial frame.
The turn rate of the sensor is given by the term . However, the gyroscopes measure , thus the term must be compensated. Particularly for pedestrian navigation, the transport rate is negligible, because the travelled distances are not large enough to significantly change its latitude and longitude within a time stamp. The Earth rotation, which is approximately 15 , is usually not compensated because the noise of the MEMS gyroscopes is much greater. Therefore, it is assumed that .
In this work, for practical reasons, we will use to shorten the term . The estimation of the Euler angles is derived from the angular rate at each time stamp .
To compute the Euler angles (see [
21]), first the biases have to be subtracted from the angular rate.
Then, the corrected turn rate measurements are integrated to compute the orientation. The orientation can be represented as a direct cosine matrix,
, which is a 3×3 rotation matrix in which each column is a unit vector along the sensor axes specified in terms of the navigation axes. The rotation matrix, the Euler angles and the quaternions are analogue ways of representing the orientation. In this work we use Euler angles instead of quaternions because the pitch angle never reaches 90
while walking. For better understanding of the range values of the pitch angle, see [
22].
The orientation at the time stamp
k is the orientation at the time stamp
rotated by the change in orientation that took place within the last
seconds, represented in a matrix form as
:
being
where
and
Therefore, the predicted Euler angles are extracted from
As Equation (
1) shows, the biases of the gyroscopes are included in the state vector together with the Euler angles. In order to present the prediction model of the biases, we first introduce our error model of the gyroscopes. The turn rate measurements,
, can be represented as
being
the error free turn rate and
the measurement error. The turn rate error can be decomposed in two errors
where
is the bias systematic error and
is the sensor noise that can be modelled as Gaussian white noise. To determine the biases error we choose an auto-regressive model of order one (AR1) [
23]. The AR1 model is defined as
The biases follow an exponentially correlated noise term defined in the AR1 model as the constant c, which is equal to the exponent , where is the correlation coefficient and can be modelled as Gaussian white noise with standard deviation . Thus, the predicted biases are
Therefore, Equation (
7) and Equation (
11) represent the vector state after the prediction stage. The values of the variance-covariance matrix of the state vector
stem from the noise of the gyroscopes through matrix
for the Euler angles and from the model presented in Equation (
10) for the biases.
2.1.2. Detectors
We want to use the accelerometer and the magnetometer to extract additional orientation information, however, these sensors should not be used continuously. We have implemented two detectors that allow using the accelerometer and the magnetometer measurements only within the right periods.
While for foot-mounted sensors the zero velocity update (ZUPT) is used [
24,
25], for pocket-mounted sensors these periods of zero velocity do not exist, thus we use the acceleration. An accelerometer is capable of measuring specific force, whose units are
. The movement of the sensor provokes accelerations that are measured as well as the gravity. In order to properly use the acceleration measurements for the orientation estimation, only the gravity acceleration has to be measured. The
zero acceleration detector identifies during the walk the periods of zero or quasi-zero acceleration due to specific movements of the sensor. Taking the advantage that the direction and magnitude of the field are known, by measuring only the gravity, it is possible to determine the roll and pitch angles of the sensor. The yaw angle cannot be determined, since it describes rotations around the axis parallel to the gravity vector [
26].
The
magnetic disturbances detector identifies during the walk the periods of constant or quasi-constant magnetic field. This is necessary because the measured magnetic field is often perturbed in indoor environments and the proposal of the state of the art we want to evaluate in this work can only be applied within these periods [
26].
2.1.3. Update Stage
There are different updates that can be applied to enhance the orientation obtained from the turn rate measurements. The authors in [
27] propose an update using the gravity field. There is also an update based on magnetic field measurements applied during constant or quasi-constant periods, which has been proposed in [
6] and further analyzed in [
16]. It is also possible to use the homogeneous magnetic field to compute the yaw angle of the sensor, but indoor and urban scenarios rarely present homogeneous magnetic fields. Turn rate measurements can also be used under the assumption that no rotation is undergoing, therefore any measured turn rate is due to gyroscopes biases. This update has been proposed in [
28]. Two updates based on acceleration and magnetic field measurements will be analyzed in detail in this work: the
Absolute Gravity Update and the
Differential Magnetic Field Update.
Absolute Gravity Update
The knowledge of the gravity field yields an estimation of the attitude if the acceleration due to the movement of the sensor is zero or quasi zero, then 9.8 . We use the zero acceleration detector to identify these periods.
Within the periods that the aforementioned condition is fulfilled, the attitude angles roll and pitch can be extracted as follows:
and
where
and
represent the Euler angles roll and pitch and
for
} represents the acceleration reading for the
i-axis measured in the sensor frame.
Then, the update equation yields:
where
is the Kalman gain [
19]. The variance-covariance matrix of
and
is equal to the mutually uncorrelated noise of the accelerometers transformed by Equations (
12) and (
13).
Differential Magnetic Field Update
This update requires a constant or quasi-constant magnetic field. We use the magnetic disturbances detector to identify these periods.
The magnetic field at the current time stamp can be computed applying the rotation of the last
seconds,
, to the magnetic field measured at the previous time stamp,
, as follows:
Then, the update equation yields:
being
the measured magnetic field and
the Kalman gain. The variance-covariance matrix of
incorporates the mutually uncorrelated noise of the magnetometers and the variance-covariance matrix of
, which contains the noise of the gyroscopes.
The
Differential Magnetic Field Update has been proposed in the state of the art in [
6,
17] to estimate the biases of the gyroscopes. The estimation of the biases is modified by this update because the biases are directly included in the update Equation (
15) through the
matrix (see Equations (
4) and (
5)). The effectiveness of this update will be evaluated in this work.
2.2. Generation of the Emulated Magnetic Field
In this section, we detail the computation needed to emulate magnetic field measurements using error-free turn rate measurements. Additionally we define the noises we add to the error-free measurements in order to obtain similar quality measurements as if they were recorded with MEMS sensors. Lastly, we validate the fiber optic gyroscopes (FOG) turn rate measurements as quasi-error-free.
We will use the Inertial Measurement Unit (IMU) DSP-1750 from KVH [
29], which embeds a FOG and a MEMS accelerometer (see
Figure 3). The acceleration measurements,
, contain noise and biases characteristic from the MEMS. The turn rate measurements of the FOG,
, are considered error-free. These error-free turn rate measurements are used to emulate different magnetic field measurements,
, as indicated in the diagram of
Figure 3.
In order to emulate MEMS gyroscopes measurements
we artificially add white Gaussian noise
and biases
to the FOG measurements as follows:
The noise and biases values have been chosen based on the typical values of the MEMS MTw from Xsens [
30] and we will detail them in the evaluation section.
To emulate magnetic measurements
the orientation of the sensor
is computed from the error-free turn rate measurements and multiplied by the magnetic field
. Finally white Gaussian noise
is added as follows:
where
is either the homogeneous magnetic field
or the perturbed magnetic field
. The
value is chosen based on the noise of the magnetometer integrated in the medium-cost MEMS MTw from Xsens and we will detail it in the evaluation section.
The election of the chosen locations on the Earth is due to their different magnetic field distributions. The homogeneous magnetic field values,
, have been rounded compared to the real values: for example, the Munich Earth Observatory has reported an average intensity of 48
T, an inclination angle of 64.23
and a declination angle of 2.57
for April 2015 [
31].
Table 1 shows all the homogeneous magnetic fields,
, expressed in East-North-Up coordinates and measured in
T.
In order to generate the perturbed magnetic field,
, we compute analytically the equations of a homogeneous field perturbed by one single ferromagnetic object. Then, we create a template of several objects. The ferromagnetic object is an iron cylinder with infinite length along the Up-axis with relative permeability
. Due to this geometry, the infinite cylinder will not influence the resulting perturbed magnetic field in the Up-axis. We generate a perturbed magnetic field from the homogeneous magnetic field of the Equator (
Table 1). The inner magnetic field of the object is not considered, since the pedestrian cannot step into the cylinder.
The magnetic field can be derived from the Maxwell’s equations, which are four partial differential equations that describe how electric and magnetic fields are generated an altered by each other and by charges and currents. If the magnetic field is static, the time derivatives are zero and the system of equations decouples. The magnetic strength vector field can be expressed as the gradient of a scalar potential field,
:
Since the influence of the ferromagnetic object on the homogeneous magnetic field decays with increasing distance, the solution of the boundary value, , at infinite distance converges to the homogeneous magnetic field , in this case the Equator.
The desired magnetic field must fulfill the boundary condition of Equation (
19), therefore the solution yields [
32]
being
the unity vectors in East and North direction respectively,
the radius of the cylinder,
the relative permeability of the material of the cylinder,
and
Using Equation (
21) and the reference magnetic field of the Equator, two different fields have been generated: one perturbed by a cylinder of
and another perturbed by a cylinder of
. We have then composed a magnetic field template with four small cylinders representing streetlights and one big cylinder representing a parked car (
Figure 4). As previously mentioned, we do not consider the influence of all cylinders simultaneously. In order to ensure that the
Differential Magnetic Field Update is correctly applied, we will not apply the update if the pedestrian walks over the transition of two fields.
FOG Turn Rate Measurements
To assess that the FOG measurements can be considered quasi-error-free, 14
of static turn rate measurements have been recorded to compute the Allan deviation shown in
Figure 5. The continuous lines depict the noise analysis of the FOG and the dashed lines depict the noise analysis of the gyroscopes of the MEMS MTw, see [
33]. The colors represent the mutually orthogonal gyroscopes in the axes
x,
y and
z.
The noise analysis of both sensors shows a decreasing trend on the left side of the plot, where the white noise is dominant. This value can be directly extracted by intersecting the curves at 1
. The FOG white noise is two orders of magnitude lower than the medium-cost MEMS gyroscopes white noise. On the right side of the plot, the Allan deviation shows a change in the trend. In this region the biases, i.e., the slow changing errors, become dominant. The bias stability,
B, is located at the minimum. The standard deviation of the biases noise is
being
where
is the averaged time corresponding to the
B and
is the sampling frequency. Therefore, the lower the
B value and the greater the
are, the more stable are the biases of the gyroscopes. The
B value is two orders of magnitude smaller for the FOG than for a medium-cost MEMS gyroscopes and
one order of magnitude greater for the FOG than for the medium-cost MEMS gyroscopes [
23]. Consequently, the FOG biases are considerably more stable than the medium-cost MEMS gyroscopes biases.
In order to evaluate the magnitude of the FOG biases and their effect on the orientation estimation, the integration of 1
of turn rate FOG measurements after correcting the Earth rotation has been computed. The FOG does not measure the transport rate since the aforementioned measurements were static, however it measures the rotation of the Earth,
, which has a value of 7.29 × 10
−5 rad s
−1, approximately 15 ° h
−1. Before integrating the turn rate measurements to compute the orientation, the Earth rotation has to be compensated as follows:
being
the transformation matrix from the navigation frame to the body frame and
where
is the latitude.
The resulting orientation angles, calculated without subtracting the biases estimation or applying updates, contain slight errors that do not exceed 1 over 1 in any angle. This error is due to the white noise and biases of the FOG which have not been compensated. A MEMS gyroscopes biases value of 0.1 ° s−1 yields to 360 over 1 . Since the FOG biases are orders of magnitude smaller and more stable than the medium-cost MEMS gyroscopes biases and the experiments last less than 1 , we consider the error due to FOG biases negligible and the FOG turn rate measurements as quasi-error-free.