Abstract
We introduce an algorithm that provides robust three-dimensional orientation of a smartphone for pedestrian indoor localization. The algorithm focuses on integration of the magnetometer and a reformulated observation model such that the influence of magnetic anomalies is mitigated. The methodological novelty of this approach lies in the use of an extended Kalman filter (EKF), based on a state vector that contains only the slow-varying systematic deviation components of the magnetometer. We apply a statistical test to the EKF residuals to detect the presence of magnetic anomalies and update the absolute heading when beneficial conditions prevail. Otherwise, the heading is propagated based on gyroscope observations. We investigate the properties of the proposed algorithm by using simulated smartphone sensor observations with different scenarios of systematic deviations. In experiments with very accurate ground truth, the proposed algorithm achieves a root mean square error of 17.4° for the computed heading, outperforming state-of-the-art algorithms by at least 40%.
1 INTRODUCTION
The computation of smartphone orientation is an important step in the process of pedestrian indoor localization, for example, when providing navigation to a certain location (Ehrlich & Blankenbach, 2018; Moder et al., 2018) or using the movement behavior of persons to gain insight into building utilization (Burgess et al., 2018; Kanda et al., 2007). We propose a new orientation-estimation algorithm (OEA) based on self-contained sensors, with a focus on magnetometer integration to provide robust absolute smartphone heading information.
The magnetometer observation model exhibits classical internal sensor errors such as those related to bias, scale factor, and misalignment (Renaudin et al., 2010) as well as platform- and environment-dependent errors such as hard-iron bias, soft-iron scaling, and magnetic anomaly bias (Groves, 2013). These latter errors are the primary reason that magnetometer integration in OEAs can cause large deviations in heading estimations. Routines for determining and mitigating these systematic deviations are mandatory, and herein, we consider only procedures that are directly applied in the indoor localization process (i.e., on-site). There are two main approaches for on-site determination or mitigation of systematic deviations, namely, instruction-based procedures and in-run procedures (Martin et al., 2016). In instruction-based procedures, the smartphone must undergo special movements or trajectories prior to the localization process; in in-run procedures, sensor errors are determined during localization. The classic instruction-based approach for smartphones is the ellipsoid-fitting approach, in which the smartphone is rotated around its three main axes to determine sensor errors. With the constraint that magnetometer observations are ideally located on a sphere with known radius, i.e., the known value of the magnitude of the earth’s magnetic field (EMF), it is possible to determine sensor and platform errors (Gebre-Egziabher et al., 2006; Klingbeil et al., 2014; Renaudin et al., 2010; Vasconcelos et al., 2011).
The obvious problem is that environment-dependent errors (i.e., magnetic anomalies) cannot be captured with the above-described approach. Thus, it is necessary to apply an in-run procedure in the OEA to determine the presence of magnetic anomalies in magnetometer observations. A common approach is to fuse the gyroscope, accelerometer, and magnetometer in an extended Kalman filter (EKF) to estimate the smartphone orientation (Gebre-Egziabher et al., 2004; Han & Wang, 2011). In these approaches, it is critical to note that magnetic anomalies also influence the inclination component of the orientation (roll φ and pitch θ in Euler angle parameterization). Madgwick et al. (2011) and Valenti et al. (2015) attempted to avoid this effect by using complementary filters, where the change in inclination is computed via the accelerometer sensor and the change in yaw (ψ in Euler angle parameterization) is computed via the magnetometer sensor. A possibility for mitigation is to use the known quantities of the EMF to determine whether magnetic anomalies are present. The literature contains many proposed solutions exploiting this information: Costanzi et al. (2016) used this knowledge to compute adaptive variances for magnetometer observations. Renaudin & Combettes (2014) and Lee et al. (2018) performed the EKF update only if conditions on the magnetometer observations were met, and Afzal et al. (2011) used statistical tests and fuzzy inference to determine the resulting heading error due to magnetic anomalies. Another approach is to parameterize the systematic deviations in an EKF (i.e., include them in the state vector). It is not common to estimate each type of systematic deviation separately in the EKF but to subsume them in a bias and an affine transformation parameterized with a symmetric matrix (Klingbeil et al., 2014; Renaudin et al., 2010). Crassidis et al. (2005) included this bias and the elements of a symmetric matrix in the state vector of an EKF and an unscented Kalman filter. Han et al. (2017) additionally integrated the EMF in the state vector and propagated it with the gyroscope observations.
Each of the previously mentioned approaches exhibits at least one of the following problems: On the one hand, systematic deviations in magnetometer observations cannot be sharply separated from the desired signal; rather, they can only be separated to a certain extent. When orientation parameters are included in the state vector of an EKF, small amounts of systematic deviations that cannot be detected influence the estimated states. On the other hand, trajectories from pedestrians do not usually contain sufficient information to decorrelate the state vector components, including orientation parameters and different types of systematic deviations.
Consequently, biased estimates arise in the presence of unmodeled effects such as magnetic anomalies. We propose an OEA that reduces these problems under some loose restrictions. The computation of the inclination, heading, and systematic deviations is divided into separate modules to avoid undesired effects due to correlated parameters. Instead of absolutely computing the heading in every epoch, the algorithm propagates the heading with the gyroscope, and the absolute part is updated periodically when good conditions for EMF sensing prevail. The update of the absolute heading is triggered by an EKF, whose state vector consists only of the magnetometer bias subsuming the systematic deviations with low temporal variance (i.e., related to the sensor and its platform). This algorithm enables the reliable detection of rapidly varying systematic deviations (i.e., magnetic anomalies and their secondary effects); therefore, the update of the absolute heading can be suppressed in such phases. However, an OEA based solely on an inertial measurement unit and magnetometer is dependent on the condition that error-free magnetometer observations are available in at least some time windows (ideally in the initial phase of the algorithm). For edge cases in which magnetometer observations are faulty all the time, it is necessary to integrate additional observation types into the OEA.
In Section 2, we introduce the proposed OEA in detail. For convenience, a list of symbols used in the extensive equations is provided in Appendix B. We validate our OEA in Section 3.1 with simulated data to show that the OEA exhibits the proposed properties. In Section 3.2, we use measured smartphone sensor data and ground truth data from a laser tracker to compare our OEA with three of the OEAs mentioned above.
2 PROPOSED ALGORITHM
The coordinate frames and smartphone sensor observations are shown in Figure 1. We assume that the observations from the gyroscope ωB, accelerometer aB, and magnetometer mB are available in the common smartphone body frame (B-frame). The aim of smartphone orientation determination is to link the B-frame with the navigation frame (N-frame), which corresponds to the coordinate system in which the pedestrian motion is described. To describe the rotations needed to link the B-frame with the N-frame, Euler angle parameterization is used herein. All frames are right-handed, and all angles are counted positive in the counter-clockwise direction. We use an intermediate frame, the local-level frame (L-frame), to separate the computation of the inclination from the computation of the heading.
Visualization of the used coordinate frames
The inclination computation (i.e., computing φ and θ) with ωB (angular velocity around the B-frame axes) and aB (components of the gravity vector gN in the B-frame) links the B-frame with the L-frame. The heading describes the rotation around the vertical z-axis and links the L-frame with the N-frame. The heading is computed by using mB (components of the EMF vector hN in the B-frame), the vertical angular velocity 
, and the magnetic declination hD. Values for hN and hD are derived from an EMF model.1
Figure 2 shows an overview of the proposed OEA for determining smartphone orientation. The inputs are the sensor observations 
, 
, and 
as well as an initial smartphone heading ψ0 at time k = 0. The algorithm consists of four modules: The computation of the inclination follows the algorithms from Särkkä et al. (2015) and Hostettler & Särkkä (2016), as outlined in Appendix A. Heading propagation is conducted with 
derived from the inclination-EKF (Section 2.1). The bias-EKF in Section 2.2 is based on a novel functional model. Its state vector contains the slow-varying systematic deviations in the magnetometer observations and uses φk, θk, ψk as well as mB as observations. In the heading update module (Section 2.3), the results from the bias-EKF are collected in time windows. If these data exhibit certain properties, the absolute heading is updated. Additionally, a one-time check of ψ0 is performed in the heading update.
Overview of the proposed OEA
It is not necessary for the user to carry the smartphone in a certain mode. The smartphone can be carried in a pocket or bag but can also be held by the user in texting or calling mode. One of the following two conditions must be met for the proposed OEA to provide useful results: Either ψ0 must be accurately known or the first couple of seconds of magnetometer observations must be free from systematic deviations (e.g., ensured by magnetometer calibration with ellipsoid fitting). We favor the second case, as it is common in several smartphone applications, requiring the user to perform certain rotations of the phone in order to trigger a built-in magnetometer calibration procedure. Still, the calibration with ellipsoid fitting must not result in an unbiased yaw angle. When ellipsoid fitting fails, there should be at least some time windows without magnetic anomalies in order to have an opportunity to improve an initial heading that is most likely erroneous. The proposed algorithm is dependent on this fact. If magnetic anomalies are present all the time, the smartphone heading can only be propagated with the relative changes derived from the gyroscope. In this case, drift effects due to the summation of gyroscpe bias must be considered, and additional observations are needed to stabilize the resulting smartphone heading. These limitations must be considered when the proposed algorithm is applied.
2.1 Heading Propagation
The heading ψ for the current epoch k is derived by integrating 
with the constant time interval between observations dt:
1
This integration is performed from the epoch κ on, at which the absolute heading has been previously updated, where dψκ|k is the accumulated heading change. The computation of ψκ is outlined in Section 2.3. The variance 
is derived with variance propagation, neglecting temporal correlations (i.e., autocorrelation). The propagated heading is the outcome of the proposed OEA but is also used as an observation in the bias-EKF, as described in the following section.
2.2 Bias-EKF
We reformulate the magnetometer model (Groves, 2013) to develop a functional model of the bias-EKF:
2
Δsn is a matrix that accounts for scaling and non-orthogonality of the magnetometer. 
and 
are the sensor bias and hard-iron bias, and 
accounts for magnetic anomalies in the sensor’s environment. Δsi contains the soft-iron effects, which are dependent on smartphone position and orientation with respect to magnetic anomalies. ∊m is the magnetometer white noise, and I is the identity matrix. 
is the rotation matrix linking the N-frame and B-frame; this term is computed from φ, θ, and ψ. We rearrange Equation (2) such that two new terms summarize the systematic deviations:
3
subsumes the deviations that are related only to the magnetometer and its carrier platform (i.e., the smartphone). This term is treated as a slow-varying bias. 
is the bias that exhibits rapid changes upon movement through the environment when magnetic anomalies are present.
As mentioned in Section 1, we choose a minimum parameterization of the state vector for the bias-EKF, which equals 
. The prediction is performed with a random walk model:
4
where 
is the predicted state, 
is the previously estimated state, and ζδ is the system noise. For better readability, we omit the index k from here on. Only the quantities related to the previously estimated state are indexed with k − 1. The filter innovation wm is computed by using mB and the Euler angles as observations:
5
Here, we use a slightly different formulation of the innovation computation that does not exhibit the common structure of “observed minus computed.” Ettlinger et al. (2018) and Vogel et al. (2018) provided a detailed explanation of this EKF variant. In the presence of magnetic anomalies, mB contains 
(see Equation (3)), which is absorbed by 
only to a certain extent depending on the preset variances of mB and ζδ (Table 1). Consequently, 
becomes biased, and the remaining influence of 
affects the residuals 
of the bias-EKF. We use 
and the corresponding variance-covariance matrix (VCM) 
to formulate a statistical test (global model test):
6
Settings of the Proposed Algorithm
σζa,i,j and 
are only relevant for the measured data containing user motion in Section 3.2.
We use the above equation to verify the compliance between the observed data and the model assumptions of the bias-EKF. In the null hypothesis H0, the expectation of the residuals is assumed to be zero 
, and therefore, TG follows a chi-square distribution with three degrees of freedom (DoFs) (i.e., the dimension of wm). In the presence of magnetic anomalies, 
also affects 
, leading to 
; therefore, TG should become significant (i.e., larger than the corresponding critical value TG,c(α), where α is the type I error). In this case (i.e., under the alternative hypothesis HA), TG follows a non-central chi-square distribution 
with non-centrality parameter λ.
Finally, we rotate mB and 
into the L-frame:
7
which we use in the following heading update step.
2.3 Heading Update
The heading update step of the proposed OEA contains two actions: the update of ψ when certain conditions are met and the one-time check of ψ0. The following equations describe the update:
8
We collect mL, 
, and TG in non-overlapping time windows with length dtW,u. To update ψκ, two conditions must be fulfilled (first line in Equation (8)). The percentage pT (second line in Equation (8)) is computed from the number of test values TG,i that exceed the corresponding critical value TG,c (α) (third line in Equation (8)). If pT > α, we assume that magnetic anomalies are present, which cannot be absorbed by 
, leading to a bias in the computed heading. The second condition dψκ|k ≤ γdψ suppresses the update of ψκ if the user is turning within a time window because, in this case, the values in the time windows are not valid for the current epoch k. γdψ is a preset threshold, and the accumulated heading change dψκ|k from Equation (1) is used as an indicator for user turns.
Up to this point, the proposed OEA relies on the the unbiasedness of ψ0, as it is used for the heading propagation and in the initialization of the bias-EKF. We use the following procedure only one time at the beginning of the trajectory to verify the compliance between ψ0 and the magnetometer observations:
9
This value is independent from the results of the bias-EKF, and if the conditions are fulfilled, ψκ is updated accordingly and the bias-EKF is re-initialized. We collect only mL in a time window with length dtW,0 starting from k = 0. First, the difference |ψk – ψm | > γψ between the propagated heading ψk from Equation (1) and ψm must exceed the predefined threshold γψ (otherwise, there is no reason for an update). The second condition is equivalent to the condition in Equation (8) and suppresses an update if the user turns within the time window. The third condition Pm > γm requires that the percentage pm (second line in Equation (9)) of 
, which fulfills two additional requirements c1 and c2 (third line in Equation (9)), is higher than the predefined threshold γm. Because we want to use the heading computed from raw magnetometer observations in this control procedure (i.e., 
is not used in the first line of Equation (9)), the observations must be checked for systematic deviations. c1 requires that the magnitude of mL be equal to the magnitude of the EMF vector hN (also derived from the EMF model in Section 2). This condition is too sparse, as all accepted solutions theoretically lie on a sphere. Thus, c2 requires that the z-component of mL be equal to the z-component of hN. All acceptable solutions of c2 lie on a plane, and the intersection with the sphere from c1 results in a circle in the three-dimensional (3D) space of solutions that fulfill both conditions. Because of measurement noise, these conditions cannot be met exactly but must be fulfilled within the threshold γσm · σm, where σm is the standard deviation of ∈m and γσm is a predefined multiplier.
In the next section, we use simulated data and experiments performed with a high-accuracy 6-DoF reference measurement system to analyze the properties of the proposed OEA. Additionally, the results from the experiments are compared with the results of three other OEAs from the literature.
3 EVALUATION
3.1 Validation with Simulated Data
In this section, we use simulated data to validate the proposed OEA in different controlled conditions, where the influence of different systematic deviations affecting the magnetometer is exactly known. The smartphone sensor observations aB, ωB, and mB are determined from a straight trajectory, as shown in Figure 3, where the orientation of the coordinate frame equals the orientation of the N-frame. The trajectory also determines the heading or yaw angle 
, and the inclination angles are constant 
, 
. Thus, the rotation matrix 
is available (~ indicates known quantities), and aB and mB are computed as follows:
10
where ∈a, ∈m, and ∈ω are the noise vectors of the accelerometer, magnetometer, and gyroscope. The components of the noise vectors are modeled as zero-mean uncorrelated Gaussian noise with a corresponding standard deviation (see Table 1). The standard deviations are derived from the smartphones used in Section 3.2 when lying static for several minutes. As we do not model stride or step accelerations, the state of the inclination-EKF only consists of 
. All settings of the proposed algorithm are summarized in Table 1.
Map with trajectories and magnetic anomalies used in the simulation scenarios
We investigate six scenarios with different systematic deviations affecting the magnetometer observations. Additionally, a baseline scenario (scenario 0) is chosen, in which no systematic deviations are introduced. Magnetic anomalies are modeled in all scenarios (except scenario 2), representing 
. The anomalies are shown in Figure 3, modeled as a magnetic dipole according to previous work (Afzal, 2011). In scenarios 1, 3, and 4, there is only one magnetic anomaly, which occurs in the middle of the trajectory, and in scenario 5, there is one anomaly at the beginning of the trajectory. Scenario 6 contains six magnetic anomalies to imitate a more realistic indoor scenario. A slow-varying bias is modeled in scenarios 2 and 3, representing 
. The bias 
is added to mB from 15.0 s on (see Figure 4(a)). A linear increase starting from 0.0 μT in the x-component is modeled from 10.0 s to 15.0 s, at which point 2.0 μT is reached. Both deviations are deterministic quantities, i.e., no noise is added. ψ0 is drawn from a Gaussian distribution with zero mean and standard deviation σψ0, according to Table 2. To analyze the influence of the accuracy of ψ0, σψ, 0 is increased in scenario 4. All scenarios are evaluated 3000 times. We investigate the results of the proposed algorithm from these six scenarios and compare them with the results for scenario 0.
Mean deviations 
for three systematic deviation scenarios with the corresponding standard deviation (dotted line) (a) Scenario 1 (b) Scenario 2 (c) Scenario 6
The orange areas in scenarios 1 and 6 indicate the influence of magnetic anomalies, the yellow area in scenario 2 indicates the influence of the constant bias, and the blue areas indicate significant global tests according to Equation (6).
Mean and Standard Deviation of 
(Deviation of Estimated 
from Reference Value) from 3000 Evaluations
The bias column indicates whether a sensor bias is added to the simulated observations, and σψ, 0 is the standard deviation of the distribution from which ψ0 is drawn.
Figure 4 shows the deviations 
of the estimated slow-varying bias from the corresponding known values for scenarios 1, 2, and 6. The results of scenario 4 do not differ from those shown in Figure 4(a), except during the first 3 s before the one-time check (Equation (9)) of ψ0 is performed. The bias-EKF exhibits the proposed behavior. 
does not absorb the magnetic anomaly, which spreads on the residuals; therefore, TG becomes significant (Figure 4(a)). An update of ψκ is not performed near the magnetic anomaly, as the heading angle would be biased and 
takes on the correct value only a couple of seconds after the magnetic anomaly is passed. In scenario 2 with the sensor bias, 
requires several seconds to follow the bias. Even if TG never becomes significant (i.e., ψκ is updated with slightly biased values in Equation (8)), 
approaches zero again after the appearance of the bias. Moreover, in the case of multiple anomalies (Figure 4(c)), the bias-EKF exhibits the desired behavior. The global test values are significant in the vicinity of magnetic anomalies and converge back to zero after approximately 5 s.
The mean deviations of the resulting heading and their standard deviations are summarized in Table 2 for all scenarios. The mean deviation is slightly increased in scenarios 2, 4, and 5 compared with the other scenarios, which provide mean deviations similar to that of scenario 0. The issue with the sensor bias in scenario 2 was discussed in the previous paragraph. Scenarios 4 and 5 represent problematic conditions at the beginning of the trajectory (i.e., poor accuracy of ψ0 or systematic deviations in mB). These two scenarios also exhibit considerably higher values of 
. This result indicates that problematic conditions in the initialization phase of the proposed OEA have the greatest impact on the achievable heading accuracy.
This simulation study is not all-encompassing, and real conditions have only been partially modeled. Yet, the results in this section indicate that the proposed OEA exhibits the desired properties. The OEA provides reasonable results for smartphone orientation when strong magnetic perturbations are present in the environment. In the next section, we analyze results from a small-scale experiment with high-accuracy ground truth values for smartphone orientation to further demonstrate the potential of the proposed OEA.
3.2 Comparison with Measured Data
We performed experiments in the measurement laboratory of the Department of Geodesy and Geoinformation, TU Wien with three smartphones (Samsung Galaxy S10, LGE Nexus 5X, and Google Pixel 5) to evaluate the achievable orientation accuracy of the proposed algorithm. The reference values of the Euler angles 
, 
, 
are determined with a Leica LTD800 laser tracker. The laser tracker provides observations in its local coordinate frame, which we denote as the Lt-frame. Pillars in the laboratory have known coordinates in a local north–east–up coordinate system, which is used to visualize the trajectories in Figure 6. The N-frame is oriented in the same way as the local north–east–up coordinate system. The coordinates of these pillars in the Lt-frame are derived from laser tracker observations to a corner cube reflector placed on the pillars. The rotation matrix 
linking the Lt-frame to the N-frame is computed with an overdetermined similarity transformation. We assembled a platform with a 3D printer that can be carried by the user and that holds the smartphone as well as the T-probe (see Figure 5). With this device, the laser tracker provides the six DoF parameters linking the T-probe coordinate frame (Tp-frame) with the Lt-frame. The T-probe is mounted with known angular offsets from the smartphone body coordinate frame (B-frame) such that the rotation matrix 
is available. The overall rotation matrix is composed from the previously described rotation steps:
11
Experimental setup
The top row of images shows the measurement area with the laser tracker and magnet. The bottom row of images shows the 3D-printed platform with the rigidly mounted T-probe and smartphone and how the platform was held during the measurements.
Trajectories 2, 3, and 4 plotted on a colormap representing the magnitude of the deviation vector ΔmB
Trajectory 1 is not plotted, as it is nearly the same as trajectory 2. The blue square represents the magnet, which is positioned in the experiment area.
The standard deviation of the angles linking the N-frame with the Lt-frame from the overdetermined similarity transformation is σφθψ N–Lt = [0.001,0.02,0.001]T [°], The standard deviation of the angles linking the Lt-frame with the Tp-frame is σφβψ, Lt – Tp = [0.002,0.002,0.002]T [°], which is determined with the T-probe being static. Because these angles are measured kinematically, we assume that σφθψ, Lt–Tp is increased by a factor of 10. Still, the uncertainty is clearly beyond a tenth of a degree for the “cumulated” angles linking the N-frame with the Tp-frame. The largest source of uncertainty is the mounting of the smartphone on the 3D-printed platform (σφθψ, Tp–B) and the realization of the B-frame in/on the smartphone. By carefully aligning the longitudinal side and back side of the smartphone with plastic screws (see Figure 5) to the 3D-printed platform, we ensure a precise alignment. From the numeric simulations in Section 3.1, we know that the accuracy of the resulting heading from the proposed algorithm is 4 − 5° under ideal conditions. We assume that this value increases by a factor of at least 2 under real conditions (platform and sensor imperfections, user motion, etc.). Thus, the chosen reference measurement system should be sufficient to provide reference values for the experiments and the following analysis.
The smartphone sensor observations aB, ωB, and mB as well as the rotation matrix 
are collected kinematically. Time synchronization between the smartphone and tracker data is conducted by using cross-correlation. The two signals for cross-correlation are ωB and 
, which is derived according to previous work (Titterton & Weston, 2004):
12
Figure 6 shows trajectories 2, 3, and 4 (trajectory 1 is the same as trajectory 2) as well as the magnet (blue square), which was positioned in the experiment area. The colormap represents the magnitude of the deviation vector 
, computed from trajectory 2 by two-dimensional interpolation. Trajectory 1 was obtained in the same manner as trajectory 2 but without the magnet. Each phone is rotated around its three main axes to trigger the hard-iron bias calibration before trajectory 1 begins. We use a magnet with a known position to have at least one magnetic anomaly influencing the magnetometer observations. As shown in Figure 6, there is another magnetic anomaly present in the upper right region of the experiment area. This pattern could be reproduced with each phone for trajectories 1, 2, and 4 (at the beginning of the trajectory).
We compute the Euler angles for all 12 trajectories with the proposed OEA. To initialize the Euler angles, we compute φ and θ with aB at time k = 0 (Appendix A), and 
is derived from Equation (11) via the laser tracker measurements. Figure 7 (top) shows the deviations of the computed EMF in the B-frame with respect to the known values for the LGE Nexus 5X and trajectory 2. Figure 7 (bottom) presents the deviations of the computed Euler angles from the corresponding known quantities. The influence of the magnet and the anomaly in the upper right region of the experiment area on the components of hB (and therefore on ψ) is very weak. Noticeable deviations of up to 25° occur in a few cases (e.g., after the initial check with Equation (9) and approximately 110 s), but they rapidly decrease to zero. The deviations of φ and θ basically have a zero mean, with a maximum of ±7°.
Deviations of hB (top) and φ, θ,ψ (bottom) for trajectory 2 and LGE Nexus 5X The durations for which the phone is closer than 1.0 m to the magnet are marked in orange (denoted “magnet” in the legend), and those for which the phone is closer than 2.0 m to the anomaly in the upper right region of the experiment area are marked in yellow (denoted “anomaly” in the legend). The turquoise line marks the initial check with Equation (9), and the purple lines mark the updates of ψκ with Equation (8).
We use three algorithms for comparison: the complementary filter from Madgwick et al. (2011) (called Madgwick or “madgw” in some plots), the EKF from Renaudin & Combettes (2014) (called Magyq), and the EKF from Han et al. (2017) (called Han). Madgwick and Magyq fuse aB, mB, and ωB in one algorithm and deliver the quaternion, which describes the 3D rotation from the N-frame into the B-frame. Thus, in contrast to the algorithm proposed herein, the computations of the inclination and heading component are not totally separated in these two approaches. Han does not estimate orientation parameters; instead, it estimates the EMF in the B-frame hB together with a magnetometer bias and symmetric matrix elements subsuming Δsn and Δsi. The estimation of the EMF (implicitly containing the heading information) and two categories of systematic deviations is the main difference of the proposed algorithm, which has minimal parameterization in the bias-EKF. The inclination angles are computed in the same way for the Han algorithm as in the proposed algorithm. We use these three algorithms to compute φ, θ, and ψ for all trajectories.
Each algorithm has at least one variable that can be used for tuning. We attempted to optimize the algorithms on the basis of the root mean square error (RMSE) of the computed heading. In the Madgwick algorithm, the maximum gyroscope measurement error can be adjusted; we set this term to 0.01° / s. Renaudin & Combettes (2014) estimated the accelerometer bias (beside the orientation quaternion q and the gyroscope quaternion bias qb,ω), which we omit herein. Thus, there are two system noise components for q and qb,ω, whose standard deviations are set to σq = 0.001 and σq,ω = 0.001 s−1. Additionally, Magyq performs static period detection and outlier rejection of aB and mB by comparing their magnitudes with the nominal values. The number of observations (aB and mB) used for static period detection is set to 20, and the corresponding thresholds are set to γ1,a = 0.3 m / s2 and γ1,m = 1.0 μT. The thresholds for outlier rejection are set to γ2,a = 1.0 m / s2 and γ2,m = 3.0 μT. Han provides one tuning variable, which is a dimensionless multiplier for controlling the influence of the system noise; this term is set to 100.0.
The resulting deviations Δψ are shown in Figure 8, and the respective RMSE values are given in Appendix B (Table B1). Madgwick performs best in trajectory 3, which exhibits the most preferable properties (no anomaly in the beginning and short duration of only 12 s). The results from the other trajectories are significantly worse. For Han, the results from trajectory 1 (without the magnet) are clearly better than those for the other three trajectories. In the provided experimental setup, heavy magnetic anomalies cannot be reasonably estimated by using functional models containing several parameters for systematic deviations (i.e., bias and symmetric matrix). Due to the very low standard deviations of the system noise obtained when using Magyq, this algorithm performs better in the short trajectories 3 and 4. Trajectory 1 contains only two very high RMSE values for the Samsung Galaxy S7 and LGE Nexus 5X. If one neglected these two values, the overall RMSE would be very close to that of the proposed OEA. The drawback of Magyq (and therefore of all algorithms that compute the 3D orientation within one EKF) is that the magnetic anomalies also influence the inclination angles (see Table 3). Because Madgwick is based on a complementary filter, the RMSE values of the inclination angles are considerably smaller, even if all smartphone observations are fused in one algorithm. The proposed OEA delivers the lowest RMSE for the computed heading given an unbiased value of ψ0 (Table 3).
Deviation of the smartphone heading computed from the known quantities Δψ for all trajectories, phones, and algorithms
RMSE Values of Euler Angles Over All Trajectories and Phones in [°]
Han and the proposed OEA use the same inclination computation procedure. Thus, these algorithms have the same RMSE values for φ and θ.
We perform the evaluations again for biased values of ψ0, as shown in Figure 9, where Δψ0 is the deviation from the correct value. The RMSE values are computed for each value Δψ0 in the same way as in Table 3. Over the whole range of Δψ0, the proposed OEA delivers the lowest RMSE values of the computed heading on average (with worse performance toward Δψ0 = −90° compared with the other algorithms). It is counterintuitive that the minima of the graphs are not exactly at the correct value of ψ0 (i.e., Δψ0 = 0°). Our explanation is that the shift in ψ0 partially compensates for systematic deviations contained in the magnetometer observations. The low variation in headings when the magnet is passed leads to additional experiment-specific effects because the summation of the magnetic anomaly and sensor-related biases is similar for all trajectories. Still, this experiment with high-accuracy ground truth values for 3D orientation reveals the potential of the proposed OEA within smartphone or pedestrian localization systems.
RMSE values of the computed heading for a biased initial heading
4 CONCLUSION
We developed an OEA that provides the absolute 3D orientation of a consumer-grade device such as a smartphone. The proposed OEA exhibits the following properties:
reliable detection of magnetic anomalies,
fast convergence of the estimated bias back to the correct value after a magnetic anomaly is passed, and
minimal/reduced influence of magnetic anomalies on the computed heading.
The limiting condition is that magnetometer observations that are free from systematic deviations must be available for at least certain time spans of the trajectory. Otherwise, the heading determination can only be performed in a relative manner by using the gyroscope or additional information from other sensors. We validated and analyzed the proposed OEA by using numerical simulations. We modeled magnetic anomalies, sensor biases, and two levels of accuracy of the initial heading in different scenarios. The scenarios with problematic conditions in the initialization phase of the proposed OEA led to slightly worse results. Nevertheless, the proposed OEA provides an estimated smartphone orientation with low deviation from the correct value.
We evaluated the performance of the proposed OEA based on experiments with high-accuracy ground truth values for 3D smartphone orientation. The proposed OEA was compared with three other algorithms from the literature that also use a magnetometer. If the initial heading is unbiased, the RMSE of the computed heading is 40% lower for the proposed OEA compared with the “second-best” algorithm over all trajectories (Tables 3 and B1). For a wide range of biased initial headings, the proposed OEA also delivers the lowest RMSE values for the computed heading (Figure 9). Although these small-scale experiments with laboratory conditions exhibit some limitations, the results indicate that the proposed OEA has potential to be included in positioning systems. In future work, the performance must beevaluated in large-scale experiments with diverse data (more DoFs in the trajectories, multiple users, more smartphone holding modes, etc.). A by-product of the proposed OEA with a potentially high benefit is a global test for magnetic anomaly detection. This approach can be used to detect and locate magnetic anomalies, which can be used as features to aid in indoor positioning, e.g., with fingerprinting.
The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Program.
HOW TO CITE THIS ARTICLE
Ettlinger, A., Wieser, A., & Neuner, H. (2024). Robust determination of smartphone heading by mitigation of magnetic anomalies. NAVIGATION, 71 (1). https://doi.org/10.33012/navi.632
Appendix A: INCLINATION-EKF
In the inclination-EKF, the gravity vector in the B-frame 
is estimated; this vector is used later to compute φ and θ. aB and ωB are smartphone sensor observations utilized as input in the inclination-EKF, in which the calibrated values delivered by the smartphone are used. ωB is the control input in the prediction of 
:
and ζg is the system noise of ḡB.
The state vector of the inclination-EKF also contains stochastic resonators (Hostettler & Särkkä, 2016; Särkkä et al., 2015), which compensate for accelerations due to user motion. We include resonators for stride 
, 
, i = x, y, z and step acceleration 
, 
, i = x, y, z for each coordinate axis (where 
are the corresponding derivatives). These terms are propagated as follows:
where f0 is the base frequency, which we set to the stride frequency, and the ζ terms are the system noise components of the stochastic resonators.
aB equals the observation vector in the filter update (i.e., measurement equation) for computing the filter innovations:
In the Euler angle representation, the roll angle φ and pitch angle θ describe the inclination; these terms are computed from 
as follows:
The angular velocity around the z-axis of the L-frame or N-frame is computed as follows:
To avoid the gimbal-lock problem 
, one can also derive the orientation quaternion from 
(Valenti et al., 2015), which describes the same rotation as φ and θ. The gimbal-lock problem does not appear in the simulation study in Section 3.1 or in the experiments in Section 3.2; thus, we utilized only the Euler angle representation.
Appendix B: ADDITIONAL MATERIAL
RMSE Values of Computed Heading in [°]
The RMSE per trajectory was computed by using the deviations of all phones and one trajectory, and the RMSE per phone was computed by using the deviations of all trajectories and one phone. The bold RMSE numbers were computed from the deviations over all phones and trajectories.
Symbols Introduced in Sections 2 and 3
Footnotes
Present address, Wiedner Hauptstraße 8-10, 1040 Vienna, Austria
↵1 World Magnetic Model (WMM), online calculator: https://www.ngdc.noaa.bov/geomag/calculators/magcalc.shtml?model=wmm, accessed: 28.05.2022
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.