Development of a Kalman Filter in the Gauss-Helmert Model for Reliability Analysis in Orientation Determination with Smartphone Sensors

Abstract

The topic of indoor positioning and indoor navigation by using observations from smartphone sensors is very challenging as the determined trajectories can be subject to significant deviations compared to the route travelled in reality. Especially the calculation of the direction of movement is the critical part of pedestrian positioning approaches such as Pedestrian Dead Reckoning (“PDR”). Due to distinct systematic effects in filtered trajectories, it can be assumed that there are systematic deviations present in the observations from smartphone sensors. This article has two aims: one is to enable the estimation of partial redundancies for each observation as well as for observation groups. Partial redundancies are a measure for the reliability indicating how well systematic deviations can be detected in single observations used in PDR. The second aim is to analyze the behavior of partial redundancy by modifying the stochastic and functional model of the Kalman filter. The equations relating the observations to the orientation are condition equations, which do not exhibit the typical structure of the Gauss-Markov model (“GMM”), wherein the observations are linear and can be formulated as functions of the states. To calculate and analyze the partial redundancy of the observations from smartphone-sensors used in PDR, the system equation and the measurement equation of a Kalman filter as well as the redundancy matrix need to be derived in the Gauss-Helmert model (“GHM”). These derivations are introduced in this article and lead to a novel Kalman filter structure based on condition equations, enabling reliability assessment of each observation.

1. Introduction

Determining the orientation in pedestrian navigation with geometric-based approaches is an essential step for positioning. There are several possibilities to mathematically parameterize the orientation. In [1] three possibilities are given: the direction-cosine-matrix, Euler angles and quaternions. Discussion on the advantages and disadvantages of these concepts can be found in [2]. In this article the Euler angles—roll, pitch and yaw—are used as it is the yaw angle that is especially necessary for calculating the 2D-Position in pedestrian navigation approaches like Pedestrian Dead Reckoning (“PDR”). The orientation can be calculated directly by the use of an accelerometer and a magnetometer respectively a gyroscope, which are nowadays integrated in most of the available smartphones. Each of these sensors is subject to specific systematic effects [3], which have to be considered in detail when integrating low-cost sensors in smartphones. These systematic deviations clearly can have a noticeable impact on the resulting position of PDR, which makes determining the orientation the critical part.

Several approaches exist to fuse and improve the calculated orientation from the individual sensors. In [4] linear combinations are used, to calculate the yaw angle of the individual results from magnetometer and gyroscope. Using a Kalman filter to fuse accelerometer, magnetometer and gyroscope measurements [5,6,7] is the most common approach. Fusing these sensors offers the opportunity to detect and correct for systematic deviations in the observations. In [8] a Kalman filter is developed which uses the observed gravity and Earth magnetic field as well as the orientation quaternion which are propagated by means of the gyroscope. If static acceleration or magnetic flux is detected, systematic sensor deviations can be determined. In [9] or [10] the gyroscope is used to detect when the user is turning. If the user walks straight, systematic sensor deviations are derived. The same approach is used in [11] to smooth the resulting trajectory. A method to minimize systematic effects in acceleration measurements resulting from the steps of pedestrians by adapting the measurement noise can be found in [2]. Using additional data can improve the ability to detect systematic effects. One possibility to retrieve additional position and attitude information is to use pictures from the smartphone camera [12]. Also building plans can be used to support the calculation of the orientation [13]. The previously mentioned approaches only use observations from one device. Especially in pedestrian navigation, much data is produced from many different smartphones, which can be used for positioning in various ways. This crowd-sourced data is mainly used to support Wi-Fi fingerprinting [14,15] by extending and updating radio maps [16,17,18]. In the context of PDR, in [19] for example, trajectories from multiple users are used to minimize the influence of magnetic perturbations inside buildings.

Combining positions from signal strength observations with PDR is done, for example, in [20]—with an adaptive Kalman filter—or in [21]—with a particle filter–to improve positioning accuracy. Such signal strength observations or “Received Signal Strength” (RSS) can be used in geometric-based approaches or in feature-based approaches and require external infrastructure. Adding such positioning techniques—for example Wi-Fi [22,23], “radio-frequency identification” (RFID) [24] or artificial magnetic fields [25]—will increase the redundancy and minimize the influence of systematic effects. Geometric- and feature-based approaches are combined in [26], where also measures are derived to quantify accuracy independently from the distribution of the observations.

This paper focuses on determining the orientation for PDR with accelerometer, magnetometer and gyroscope in a Kalman filter. The focus is on analyzing the reliability of the observations from smartphone-sensors, which are non-linearly included in condition equations. The Kalman filter is commonly used to estimate states in a Gauss-Markov model (“GMM”—model with observation equations), wherein the observations of the measurement equation are linear and can be formulated as functions of the states [27,28]. However, the equations for the orientation determination in PDR are condition equations, containing implicit relations between states and observations. Thus, reliability analysis based on condition equations, necessitates a new derivation of the system equation and the measurement equation of a Kalman filter as well as the redundancy matrix in the Gauss-Helmert model (“GHM”—model with condition equations). The formulation of the system equation in the GHM can be found in [29,30], wherein the measurement equation and corresponding redundancy matrix are still assumed to satisfy the GMM. Hence, an important novelty of this article is the derivation of the whole Kalman filter structure in the GHM, enabling the possibility to calculate reliability measures for observations which are non-linearly included in the condition equations.

Section 2 contains the explanation of the existing approach and subsequent problem description of reliability, as well as the derivation of the update-equations and redundancy matrix of the reformulated Kalman filter. The data of a measured trajectory will be used to estimate the Euler angles as well as the reliability measures. These results are presented in Section 3 whereas Section 4 concludes this contribution.

2. Orientation Determination

2.1. Existing Approach

In the Kalman filter used for orientation determination in this article, the state parameters are the Euler angles roll ϕ, pitch θ and yaw ψ. Based on the assumption that roll and pitch are constant while the user is walking, they are predicted at actual epoch k with a random-walk model (1) and (2). The yaw angle is predicted with two different system equations, whereby (3) will be used if the user is walking straight and (4) will be used at detected turns. In the second case the observed turn rates from gyroscope g_y,k and g_z,k are included as control quantities summed in the vector u_k. By doing this, the predicted yaw should immediately follow the user’s turn. The noise components of the state parameters are c_k−1,ϕ, c_k−1,θ, c_k−1,ψ and Δt is the time interval between two consecutive Kalman updates:

$$f_1:{\phi }_k={\phi }_{k-1}+c_{k-1,\phi }$$

(1)

$$f_2:{\theta }_k={\theta }_{k-1}+c_{k-1,\theta }$$

(2)

$$f_3:{\psi }_k={\psi }_{k-1}+c_{k-1,\psi }$$

(3)

$$f_3:{\psi }_k={\psi }_{k-1}+\frac{\Delta t}{\cos {\theta }_{k-1}}\left(g_{y,k}\sin {\phi }_{k-1}+g_{z,k}\cos {\phi }_{k-1}\right)+c_{k-1,\psi }$$

(4)

Turn detection is done by applying a statistical test on the filter innovations of the yaw angle d_ψ,i [29,31]. The null hypothesis H₀ of this statistical test says that the vector d_ψ containing the filter innovations d_ψ,i from the last n epochs is equal to the zero vector (5), whereas the alternate hypothesis H_A states that d_ψ is significantly different to the zero vector (6). If the test value exceeds the corresponding quantile of the chi-square distribution (7), (4) will be used to predict the yaw angle. D_ψ is the variance-covariance matrix (VCM) of the innovation vector d_ψ and only contains variances on the diagonal belonging to the corresponding d_ψ,i. This is a simplification, as auto-covariances may be present. The use of random-walk (3) results in smoother trajectories in sections when the user walks straight. By neglecting the observations from gyroscope in the random-walk model, the influence of systematic sensor deviations (gyro-drift) on the filter result is minimized:

$$H_0:E({\mathit{d}}_{\psi })=E\left({\left[\begin{array}{cccc}{d}_{\psi ,1}& d_{\psi ,2}& \cdots & d_{\psi ,n}\end{array}\right]}^T\right)=\mathbf{0}$$

(5)

$$H_A:E({\mathit{d}}_{\psi })\ne \mathbf{0}$$

(6)

$$P\{{\mathit{d}}_{\psi }^T{\mathit{D}}_{\psi }^{-1}{\mathit{d}}_{\psi }\le {\chi }_{n;1-\alpha }^2|H_0\}=1-\alpha$$

(7)

For the update equations also the VCM of the predicted state is needed (8):

$${\sum }_{\overline{x}\overline{x},k}={\mathit{T}}_{k,k-1}{\sum }_{\widehat{x}\widehat{x},k-1}{\mathit{T}}_{k,k-1}^T+{\mathit{U}}_{k,k-1}{\sum }_{uu,k}{\mathit{U}}_{k,k-1}^T+{\mathit{C}}_{k,k-1}{\sum }_{cc,k}{\mathit{C}}_{k,k-1}^T$$

(8)

Therein, ${\Sigma }_{ii,k}$—with the index i = ${\widehat{x}}_{k-1}$, u, c—are the VCMs of the corresponding observation groups in epoch k. T_k,k−1 is the state transition matrix, U_k,k−1 the control matrix and C_k,k−1 is the noise matrix, each referred to epoch k. These system matrices are Jacobi matrices and in general contain the derivatives of the system equations with respect to the corresponding observation group [27,28,29]. Equations (9)–(11) show the system matrices for the approach presented in this section where E is the identity matrix.

$${\mathit{T}}_{k,k-1}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ \frac{\partial f_3}{\partial {\phi }_k}& \frac{\partial f_3}{\partial {\theta }_k}& 1\end{array}\right]$$

(9)

$${\mathit{U}}_{k,k-1}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& \frac{\partial f_3}{\partial g_{y,k}}& \frac{\partial f_3}{\partial g_{z,k}}\end{array}\right]$$

(10)

$${\mathit{C}}_{k,k-1}=\Delta t\cdot {\mathit{E}}_{3\times 3}$$

(11)

The observations in the measurement equation are directly observed Euler angles, as also shown in [32]. Therefore, the design matrix A_m,k of the Kalman filter equals the identity matrix. They are calculated outside of the Kalman filter ((12)–(14)) by using observed accelerations a_x,k, a_y,k, a_z,k and magnetic flux densities m_x,k, m_y,k, m_z,k [33]. The accelerations are filtered in a separate Kalman filter to remove high-frequency components due to the movement of the user [2]:

$$f_4:{\phi }_k={\tan }^{-1}\left(\frac{a_{y,k}}{a_{z,k}}\right)$$

(12)

$$f_5:{\theta }_k={\tan }^{-1}\left(\frac{-a_{x,k}}{a_{y,k}\sin \left({\phi }_k\right)+a_{z,k}\cos \left({\phi }_k\right)}\right)$$

(13)

$$f_6:{\psi }_k={\tan }^{-1}\left(\frac{m_{z,k}\sin \left({\phi }_k\right)-m_{y,k}\cos \left({\phi }_k\right)}{m_{x,k}\cos \left({\theta }_k\right)+m_{y,k}\sin \left({\theta }_k\right)\sin \left({\phi }_k\right)+m_{z,k}\sin \left({\theta }_k\right)\cos \left({\phi }_k\right)}\right)$$

(14)

The yaw angle calculated with (14) can be subject to magnetic perturbations, especially inside buildings. The redundant determination of the Euler angles in the Kalman filter by using the system equation and the measurement equation dampens the influence of magnetic perturbations. Additionally, the standard deviation of the magnetometer measurements σ_m will be increased if the magnitude of the measured magnetic field $\Vert {\mathit{m}}_k\Vert$ is not stable (15), leading to an adaptive standard deviation σ_m,k for each Kalman filter epoch. As the geomagnetic field should be constant related to the dimensions of a building, it will be assumed that magnetic perturbations are present if the measured magnitude changes:

$${\sigma }_{m,k}={\sigma }_m+|\Vert {\mathit{m}}_k\Vert -\Vert {\mathit{m}}_{k-1}\Vert |$$

(15)

Using the covariance propagation law, the VCM belonging to the directly observed Euler angles ${\Sigma }_{\phi \theta \psi ,k}$ is derived (16), wherein H_k is a Jacobi matrix containing the derivatives of (12)–(14) with respect to the accelerometer and magnetometer measurements.

$$\begin{array}{l}{\sum }_{\phi \theta \psi ,k}={\mathit{H}}_k\left[\begin{array}{cc}{\sum }_{aa,k}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& {\sigma }_{m,k}{\mathit{E}}_{3\times 3}\end{array}\right]{\mathit{H}}_k^T\\ {\mathit{H}}_k=\left[\begin{array}{cccccc}0& \frac{\partial f_4}{\partial a_{y,k}}& \frac{\partial f_4}{\partial a_{z,k}}& 0& 0& 0\\ \frac{\partial f_5}{\partial a_{x,k}}& \frac{\partial f_5}{\partial a_{y,k}}& \frac{\partial f_5}{\partial a_{z,k}}& 0& 0& 0\\ \frac{\partial f_5}{\partial a_{x,k}}& \frac{\partial f_5}{\partial a_{y,k}}& \frac{\partial f_5}{\partial a_{z,k}}& \frac{\partial f_5}{\partial m_{x,k}}& \frac{\partial f_5}{\partial m_{y,k}}& \frac{\partial f_5}{\partial m_{z,k}}\end{array}\right]\end{array}$$

(16)

2.2. Problem Description

To determine the user’s position, the step length is also estimated and will be used with fixed variance in PDR. Figure 1 shows a trajectory calculated with PDR, whereby the yaw is estimated with the Kalman filter mentioned in Section 2.1. Additionally, the 95% confidence ellipses as well as the reference trajectory are part of Figure 1. This measured trajectory will also be used to analyze the partial redundancy in Section 3. During the measurements, the user held the smartphone in portrait mode (Figure 2 left, ϕ_k ~ 0°). The reference trajectory is calculated from the measurements of the TS16 total station (Leica Geosystems, Heerbrugg, Switzerland) which is tracking the user by the help of a 360°-mini-prism on a helmet (Figure 2 right).

In PDR, the yaw angle is responsible for the shape of the trajectory and the step length for the scale. Hence in Figure 1, the estimated orientation causes the deviations between the estimated and the reference trajectory. As these deviations are not captured by the confidence ellipses—which are a measure for the precision [34]—there are two possible reasons for their appearance: systematic deviations in the observed data or non-Gaussian distributed data. Because of the obvious systematics in the estimated trajectory, it is assumed that systematic deviations are present in the observed data.

Reliability theory deals with the detection of large systematic deviations (inner reliability) and their effect on the estimated quantities (outer reliability) [35,36]. To identify the measurements responsible for the systematic deviations of the orientation between estimated and reference trajectory, the inner reliability is used. Here the partial redundancies r_i play a key role. According to [29,37,38], all quantities with stochastic information can be treated as observations.

Observations related to the system equation are the estimated state from the previous epoch ${\widehat{\mathit{x}}}_{k-1}$, the control variables u and the noise variables c. The observations of the measurement equation are summed in the vector ${\mathit{l}}_m$, where m labels quantities of the measurement equation. With (17)–(20) the partial redundancy can be calculated for the previously mentioned observations according to [29]:

$${\displaystyle \sum _{j=1}^{n_x}{r}_{{\widehat{x}}_{k-1},k,j}=}{\displaystyle \sum _{j=1}^{n_x}{\mathit{e}}_j^T{\sum }_{\widehat{x}\widehat{x},k-1}{\mathit{T}}_{k,k-1}^T{\mathit{A}}_{m,k}^T{\mathit{D}}_k^{-1}{\mathit{A}}_{m,k}{\mathit{T}}_{k,k-1}{\mathit{e}}_j}$$

(17)

$${\displaystyle \sum _{j=1}^{n_u}{r}_{u,k,j}=}{\displaystyle \sum _{j=1}^{n_u}{\mathit{e}}_j^T{\sum }_{uu,k}{\mathit{U}}_{k,k-1}^T{\mathit{A}}_{m,k}^T{\mathit{D}}_k^{-1}{\mathit{A}}_{m,k}{\mathit{U}}_{k,k-1}{\mathit{e}}_j}$$

(18)

$${\displaystyle \sum _{j=1}^{n_w}{r}_{c,k,j}=}{\displaystyle \sum _{j=1}^{n_w}{\mathit{e}}_j^T{\sum }_{cc,k}{\mathit{C}}_{k,k-1}^T{\mathit{A}}_{m,k}^T{\mathit{D}}_k^{-1}{\mathit{A}}_{m,k}{\mathit{C}}_{k,k-1}{\mathit{e}}_j}$$

(19)

$${\displaystyle \sum _{j=1}^{n_l}{r}_{l_m,k,j}=}{\displaystyle \sum _{j=1}^{n_l}{\mathit{e}}_j^T{\sum }_{ll,m,k}{\mathit{D}}_k^{-1}{\mathit{e}}_j}$$

(20)

Therein, n_i is the number of observations in each observation group and D_k is the VCM of the filter innovation [27,28,29]. e_j is the unity vector to select the corresponding r_i respectively diagonal element. In the above formulation of the Kalman filter, the r_i can only be calculated for the directly observed angles in the measurement equation but not for the original observations from accelerometer and magnetometer due to the structure of (12)–(14).

2.3. Kalman Filter in the Gauss-Helmert Model

In the GMM, the true observations $\tilde{\mathit{l}}$ can be modeled as a function of the true parameters $\tilde{\mathit{x}}$—which also holds for the estimated observations $\widehat{\mathit{l}}$ and estimated parameters $\widehat{\mathit{x}}$ (21). As mentioned above, the measurement Equations (12)–(14) have another structure, which matches the GHM [39] shown in (22). Hence, the update equations of the Kalman filter have to be derived in the GHM, to directly use the observations from smartphone-sensors. Afterwards, the corresponding r_i need to be derived for this case. In [30] the GHM is also used to estimate variance components for the system noise. Though, there is still the assumption of using the GMM in the measurement equation and the results cannot be used in this article.

$$\tilde{\mathit{l}}-\mathit{f}\left(\tilde{\mathit{x}}\right)=\widehat{\mathit{l}}-\mathit{f}\left(\widehat{\mathit{x}}\right)=\mathbf{0}$$

(21)

$$\mathit{f}\left(\tilde{\mathit{l}},\tilde{\mathit{x}}\right)=\mathit{f}\left(\widehat{\mathit{l}},\widehat{\mathit{x}}\right)=\mathbf{0}$$

(22)

As the functional relations (22) are non-linear in general—and especially in this article—they have to be linearized. Neither the real observations and parameters, nor their estimated values are known a priori. Hence, to do Taylor linearization, approximate values ${\mathit{l}}_0$ and ${\mathit{x}}_0$ have to be used which also have to satisfy the non-linear relations [34]. In least-squares, especially ${\mathit{x}}_0$ is assumed to be non-stochastic. According to [40], ${\mathit{l}}_0$ can be replaced by the observed data $\mathit{l}$ after linearization, which results in the linearized, functional model (23):

$$\mathit{f}(\mathit{l},{\mathit{x}}_0)+\mathit{A}(\widehat{\mathit{x}}-{\mathit{x}}_0)+\mathit{B}(\widehat{\mathit{l}}-\mathit{l})=\mathbf{0}$$

(23)

The first term corresponds to the misclosure vector w, the bracketed expression in the second term equals the stochastic additions $\Delta \widehat{\mathit{x}}$ to the approximate parameters ${\mathit{x}}_0$ and the bracketed expression in the third term equals the residuals $\mathbf{v}$. Formula (23) is only valid for the first iteration of least-squares estimation, because linearization of the functional model necessitates an iterative approach. In the subsequent iterations, the functional model is always linearized at the previously estimated observations and parameters [40]. Formulas (24)–(27) show, how the searched quantities with corresponding VCM are calculated in the GHM in the context of least-squares [41]:

$${\sum }_{\Delta \widehat{x}\Delta \widehat{x}}={\left({\mathit{A}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\mathit{A}\right)}^{-1}$$

(24)

$$\Delta \widehat{\mathit{x}}=-{\sum }_{\Delta \widehat{x}\Delta \widehat{x}}{\mathit{A}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\mathit{w}$$

(25)

$$\mathit{v}=-{\sum }_{ll}{\mathit{B}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\left(\mathit{E}-\mathit{A}{\sum }_{\Delta \widehat{x}\Delta \widehat{x}}{\mathit{A}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\right)\mathit{w}$$

(26)

$${\sum }_{vv}={\sum }_{ll}{\mathit{B}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\left(\mathit{E}-\mathit{A}{\sum }_{\Delta \widehat{x}\Delta \widehat{x}}{\mathit{A}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\right)\mathit{B}{\sum }_{ll}$$

(27)

B is the observation matrix, ${\Sigma }_{\Delta \widehat{x}\Delta \widehat{x}}$ is the VCM of the estimated additions, ${\Sigma }_{ll}$ the VCM of the observations and ${\Sigma }_{\mathit{vv}}$ is the VCM of the residuals. It is important that A and B are column-regular matrices. Otherwise the inverse matrices cannot be calculated. To avoid singular matrices, the functional relations have to be chosen, such that there will be no linearly dependent columns in these matrices.

The Kalman filter is based on sequential least-squares with an additional state transition model. The state transition model respectively system equation is described in general by non-linear, stochastic, vector-matrix differential equations. We assume, that such differential equations have the form shown in (28) after linearization, which is also basis for deriving the system equation in [27,28,29]:

$$\dot{\mathit{x}}(t)=\mathit{F}\mathit{x}(t)+\mathit{L}\mathit{u}(t)+\mathit{G}\mathit{c}(t)$$

(28)

t is the continuous time variable. The system matrix F, control-input matrix L and noise-input matrix G are assumed to be time-invariant. Solving such differential equations in the state space—also shown in [27,28,29]—unambiguously defines the state parameters at time k and gives the approximate formulation (29):

$${\mathit{x}}_k\approx {\mathit{T}}_{k,k-1}{\mathit{x}}_k+{\mathit{U}}_{k,k-1}{\mathit{u}}_k+{\mathit{C}}_{k,k-1}{\mathit{c}}_k$$

(29)

As mentioned in [27], the predicted parameter vector ${\overline{\mathit{x}}}_k$ should be calculated from the original set of functions (30) in the non-linear case to avoid linearization errors. Equation (29) is necessary for calculating the VCM of ${\overline{\mathit{x}}}_k$ with the covariance propagation law (8).

$${\overline{\mathit{x}}}_k=\mathit{f}({\widehat{\mathit{x}}}_{k-1},{\mathit{u}}_k,{\mathit{c}}_k)$$

(30)

The system can also be described by several sub-systems which yields in an over-determined system equation. This will be dealt with in the next section, where one yaw angle will be estimated from multiple trajectory data. Thus, the system equation decomposes into two sets of equations. The first set unambiguously defines the predicted parameters by using ${\widehat{\mathit{x}}}_{k-1}$, u_k and c_k which is already shown in (29) respectively (30). This set of equations will be—analogue to [29]—formulated in the GHM, such that it contains the residuals belonging to ${\widehat{\mathit{x}}}_{k-1}$, u_k and c_k (31). The derivation of (31) can be found in the Appendix A:

$${\overline{\mathit{x}}}_k-{\widehat{\mathit{x}}}_k+{\mathit{v}}_{\overline{x},k}=-\mathit{E}\Delta {\widehat{\mathit{x}}}_k+\left[\begin{array}{ccc}{\mathit{T}}_{k,k-1}& {\mathit{U}}_{k,k-1}& {\mathit{C}}_{k,k-1}\end{array}\right]\left[\begin{array}{c}{\mathit{v}}_{{\widehat{x}}_{k-1},k}\\ {\mathit{v}}_{u,k}\\ {\mathit{v}}_{c,k}\end{array}\right]=\mathbf{0}$$

(31)

The second set consists of condition equations, having the same structure like (23). These condition equations contain the same quantities like the first set of equations and can therefore be formulated as shown in (32). w^*_s,k are the misclosures, arising because of the overdetermined system equation respectively condition equations. A^*_s,k is the design matrix and T^*_k.k−1, U^*_k,k−1 and C^*_k,k−1 are the observation matrices which belong to the condition equations of the system equation:

$${\mathit{w}}_{s,k}^{*}+{\mathit{A}}_{s,k}^{*}\Delta {\widehat{\mathit{x}}}_k+\left[\begin{array}{ccc}{\mathit{T}}_{k,k-1}^{*}& {\mathit{U}}_{k,k-1}^{*}& {\mathit{C}}_{k,k-1}^{*}\end{array}\right]\left[\begin{array}{c}{\mathit{v}}_{{\widehat{x}}_{k-1},k}\\ {\mathit{v}}_{u,k}\\ {\mathit{v}}_{c,k}\end{array}\right]=\mathbf{0}$$

(32)

The functional model of the measurement equation has the structure (23). Fusing the system equation and the measurement equation leads to the functional model of a Kalman filter formulated in the GHM with additional condition equations in the system Equation (33):

$$\left[\begin{array}{c}\begin{array}{c}\mathbf{0}\\ {\mathit{w}}_{s,k}^{*}\end{array}\\ {\mathit{w}}_{m,k}\end{array}\right]+\left[\begin{array}{c}\begin{array}{c}-\mathit{E}\\ {\mathit{A}}_{s,k}^{*}\end{array}\\ {\mathit{A}}_{m,k}\end{array}\right]\Delta {\widehat{\mathit{x}}}_k+\left[\begin{array}{c}{\mathit{T}}_{k,k-1}\\ {\mathit{T}}_{k,k-1}^{*}\\ \mathbf{0}\end{array}\begin{array}{c}{\mathit{U}}_{k,k-1}\\ {\mathit{U}}_{k,k-1}^{*}\\ \mathbf{0}\end{array}\begin{array}{c}{\mathit{C}}_{k,k-1}\\ {\mathit{C}}_{k,k-1}^{*}\\ \mathbf{0}\end{array}\begin{array}{c}\mathbf{0}\\ \mathbf{0}\\ {\mathit{B}}_{m,k}\end{array}\right]\left[\begin{array}{c}{\mathit{v}}_{{\widehat{x}}_{k-1},k}\\ {\mathit{v}}_{u,k}\\ {\mathit{v}}_{c,k}\\ {\mathit{v}}_{l_m,k}\end{array}\right]=\underset{\mathit{w}}{\underbrace{\left[\begin{array}{c}{\mathit{w}}_{s,k}\\ {\mathit{w}}_{m,k}\end{array}\right]}}+\underset{\mathit{A}}{\underbrace{\left[\begin{array}{c}{\mathit{A}}_{s,k}\\ {\mathit{A}}_{m,k}\end{array}\right]}}\Delta {\widehat{\mathit{x}}}_k+\underset{\mathit{B}}{\underbrace{\left[\begin{array}{cc}{\mathit{B}}_{s,k,k-1}& \mathbf{0}\\ \mathbf{0}& {\mathit{B}}_{m,k}\end{array}\right]}}\underset{\mathit{v}}{\underbrace{\left[\begin{array}{c}{\mathit{v}}_{l_s,k}\\ {\mathit{v}}_{l_m,k}\end{array}\right]}}=\mathbf{0}$$

(33)

Now the formulas of the GHM in least-squares can be applied and the parameters with corresponding VCM can be calculated. Equation (34) shows the VCM of the estimated parameters, which is derived by inserting A and B into (24):

$${\sum }_{\widehat{x}\widehat{x},k}={\left({\mathit{A}}_{s,k}^T{\left({\mathit{B}}_{s,k,k-1}{\sum }_{ll,s,k}{\mathit{B}}_{s,k,k-1}^T\right)}^{-1}{\mathit{A}}_{s,k}+{\mathit{A}}_{m,k}^T{\left({\mathit{B}}_{m,k}{\sum }_{ll,m,k}{\mathit{B}}_{m,k}^T\right)}^{-1}{\mathit{A}}_{m,k}\right)}^{-1}={\left({\sum }_{\widehat{x}\widehat{x},s,k}^{-1}+{\sum }_{\widehat{x}\widehat{x},m,k}^{-1}\right)}^{-1}$$

(34)

${\Sigma }_{\widehat{x}\widehat{x},s,k}$ is the VCM of the parameters which would be the result of a Kalman filter only considering the system equation. Similarly, ${\Sigma }_{\widehat{x}\widehat{x},m,k}$ is the VCM of the parameters if only the measurement equation would be considered in the Kalman filter. ${\Sigma }_{\widehat{x}\widehat{x},k}$ is the inverse of the sum of these two matrices. To further process (34), the Woodbury formula for matrix inversion—according to [42]—is applied (35), where M, N, O and P are arbitrary matrices and not related to the derivations made in this section:

$${\left(\mathit{M}+\mathit{N}\mathit{O}\mathit{P}\right)}^{-1}={\mathit{M}}^{-1}-{\mathit{M}}^{-1}\mathit{N}{\left({\mathit{O}}^{-1}+\mathit{P}{\mathit{M}}^{-1}\mathit{N}\right)}^{-1}\mathit{P}{\mathit{M}}^{-1}$$

(35)

Depending on which term of the sum in (34) is chosen to be the matrix M in the Woodbury formula (35), results will show two equivalent representations (36) and (37) of the VCM of the estimated parameters. It has to be mentioned, that this VCM corresponds to ${\widehat{\mathit{x}}}_k$ and not to $\Delta {\widehat{\mathit{x}}}_k$. The reason is that ${\mathit{x}}_0$ equals ${\overline{\mathit{x}}}_k$, which is stochastic and its stochastic information is implicitly integrated by adding its calculation to the functional model ((31) and (33)). Whereas in least-squares ${\mathit{x}}_0$ is non-stochastic and therefore ${\Sigma }_{\widehat{x}\widehat{x}}$ equals ${\Sigma }_{\Delta \widehat{x}\Delta \widehat{x}}$.

$${\sum }_{\widehat{x}\widehat{x},k}=\left(\mathit{E}-{\sum }_{\widehat{x}\widehat{x},s,k}{\mathit{A}}_{m,k}^T{\left({\mathit{B}}_{m,k}{\sum }_{ll,m,k}{\mathit{B}}_{m,k}^T+{\mathit{A}}_{m,k}{\sum }_{\widehat{x}\widehat{x},s,k}{\mathit{A}}_{m,k}^T\right)}^{-1}{\mathit{A}}_{m,k}\right){\sum }_{\widehat{x}\widehat{x},s,k}$$

(36)

$${\sum }_{\widehat{x}\widehat{x},k}=\left(\mathit{E}-{\sum }_{\widehat{x}\widehat{x},m,k}{\mathit{A}}_{s,k}^T{\left({\mathit{B}}_{s,k,k-1}{\sum }_{ll,s,k}{\mathit{B}}_{s,k,k-1}^T+{\mathit{A}}_{s,k}{\sum }_{\widehat{x}\widehat{x},m,k}{\mathit{A}}_{s,k}^T\right)}^{-1}{\mathit{A}}_{s,k}\right){\sum }_{\widehat{x}\widehat{x},m,k}$$

(37)

By comparing (36) with the update equation for the VCM of the estimated parameters in the GMM (see [29] or [27]), the VCM of the filter innovation D_m,k and the gain matrix K_m,k—belonging to the measurement equation—can be found ((39) and (41)). In the same manner, the VCM of the filter innovation D_s,k of and the gain matrix K_s,k belonging to the system equation are derived ((38) and (40)):

$${\mathit{D}}_{s,k}={\mathit{B}}_{s,k,k-1}{\sum }_{ll,s,k}{\mathit{B}}_{s,k,k-1}^T+{\mathit{A}}_{s,k}{\sum }_{\widehat{x}\widehat{x},m,k}{\mathit{A}}_{s,k}^T$$

(38)

$${\mathit{D}}_{m,k}={\mathit{B}}_{m,k}{\sum }_{ll,m,k}{\mathit{B}}_{m,k}^T+{\mathit{A}}_{m,k}{\sum }_{\widehat{x}\widehat{x},s,k}{\mathit{A}}_{m,k}^T$$

(39)

$${\mathit{K}}_{s,k}={\sum }_{\widehat{x}\widehat{x},m,k}{\mathit{A}}_{s,k}^T{\mathit{D}}_{s,k}^{-1}$$

(40)

$${\mathit{K}}_{m,k}={\sum }_{\widehat{x}\widehat{x},s,k}{\mathit{A}}_{m,k}^T{\mathit{D}}_{m,k}^{-1}$$

(41)

The estimated additions to the approximate parameters are derived by using (25) and (33). Formula (42) shows these additions, where the results from the previous paragraph are already considered. Hence, it includes corrections respectively updates for the system equation as well as for the measurement equation:

$$\Delta {\widehat{\mathit{x}}}_k=-{\mathit{K}}_{s,k}{\mathit{w}}_{s,k}-{\mathit{K}}_{m,k}{\mathit{w}}_{m,k}=-\left[\begin{array}{cc}{\mathit{K}}_{s,k}& {\mathit{K}}_{m,k}\end{array}\right]\left[\begin{array}{c}{\mathit{w}}_{s,k}\\ {\mathit{w}}_{m,k}\end{array}\right]$$

(42)

If the system equation is not overdetermined, (32) will not be used in the functional model (33), which leads to simplifications in the update equations. The classical Kalman filter update equations as well as the two GHM variants are summarized in

Table 1. Comparison of the Kalman filter update equations.

	Gauss-Markov	Gauss-Helmert (Simplification)	Gauss-Helmert
$\mathit{d}$	${\mathit{l}}_m-{\mathit{f}}_m\left(\overline{\mathit{x}}\right)$	$-{\mathit{w}}_m=-{\mathit{f}}_m\left({\mathit{l}}_m,\overline{\mathit{x}}\right)$	$\left[\begin{array}{c}-{\mathit{w}}_s\\ -{\mathit{w}}_m\end{array}\right]=\left[\begin{array}{c}-{\mathit{f}}_s\left({\mathit{l}}_s,\overline{\mathit{x}}\right)\\ -{\mathit{f}}_m\left({\mathit{l}}_m,\overline{\mathit{x}}\right)\end{array}\right]$
$\mathit{D}$	${\mathit{A}}_m{\sum }_{\overline{x}\overline{x}}{\mathit{A}}_m^T$	${\mathit{B}}_m{\sum }_{ll,m}{\mathit{B}}_m^T+{\mathit{A}}_m{\sum }_{\overline{x}\overline{x}}{\mathit{A}}_m^T$	$\left[\begin{array}{cc}{\mathit{D}}_s& \mathbf{0}\\ \mathbf{0}& {\mathit{D}}_m\end{array}\right]=$ $=\left[\begin{array}{cc}{\mathit{B}}_s{\sum }_{ll,s}{\mathit{B}}_s^T+{\mathit{A}}_s{\sum }_{\widehat{x}\widehat{x},m}{\mathit{A}}_s^T& \mathbf{0}\\ \mathbf{0}& {\mathit{B}}_m{\sum }_{ll,m}{\mathit{B}}_m^T+{\mathit{A}}_m{\sum }_{\widehat{x}\widehat{x},s}{\mathit{A}}_m^T\end{array}\right]$
$\mathit{K}$	${\sum }_{\overline{x}\overline{x}}{\mathit{A}}_m^T{\mathit{D}}^{-1}$		$\left[\begin{array}{cc}{\mathit{K}}_s& {\mathit{K}}_m\end{array}\right]=\left[\begin{array}{cc}{\sum }_{\widehat{x}\widehat{x},m}{\mathit{A}}_s^T{\mathit{D}}_s^{-1}& {\sum }_{\widehat{x}\widehat{x},s}{\mathit{A}}_m^T{\mathit{D}}_m^{-1}\end{array}\right]$
${\sum }_{\widehat{x}\widehat{x}}$	$\left(\mathit{E}-\mathit{K}{\mathit{A}}_m\right){\sum }_{\overline{x}\overline{x}}$		$\left(\mathit{E}-{\mathit{K}}_m{\mathit{A}}_m\right){\sum }_{\widehat{x}\widehat{x},s}=\left(\mathit{E}-{\mathit{K}}_s{\mathit{A}}_s\right){\sum }_{\widehat{x}\widehat{x},m}$
$\widehat{\mathit{x}}$	$\overline{\mathit{x}}+\mathit{K}\mathit{d}$

.

Approaches whose aim is to detect systematic deviations in the observations respectively to quantify the inner reliability in the GMM, are often based on the disturbed residuals $\overline{\mathit{v}}$ [43]. The error-term $\nabla \mathbf{v}$ in (43) is caused by the systematic deviations in the observations, summed in the error-vector $\nabla \mathit{l}$:

$$\overline{\mathit{v}}=\mathit{v}+\nabla \mathit{v}$$

(43)

r_i is the factor which specifies how an observation deviation $\nabla l_i$ influences the corresponding residual ${\mathit{v}}_i$ (44). Hence, high r_i are desirable to detect systematic deviations [29]:

$$\nabla v_i=-r_i\nabla l_i$$

(44)

Transferring these thoughts to the GHM, (26) has to be used to derive R, containing the r_i on its diagonal. R is an idempotent matrix, whose trace equals the overall redundancy of the estimation problem [34]. Formula (26) only contains the misclosures w, which can be linearized by Bl according to [34]. If $\nabla \mathit{l}$ is taken into account, results show the disturbed model (45):

$$\overline{\mathit{v}}=-{\sum }_{ll}{\mathit{B}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\left(\mathit{E}-\mathit{A}{\sum }_{\Delta \widehat{x}\Delta \widehat{x}}{\mathit{A}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\right)\mathit{B}\left(\mathit{l}+\nabla \mathit{l}\right)$$

(45)

Thus, there is a direct relation between observations and residuals and the redundancy matrix for the GHM is now available (46):

$$\mathit{R}={\sum }_{ll}{\mathit{B}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\left(\mathit{E}-\mathit{A}{\sum }_{\widehat{x}\widehat{x}}{\mathit{A}}^T{\left(\mathit{B}{\sum }_{ll}{\mathit{B}}^T\right)}^{-1}\right)\mathit{B}$$

(46)

By using (27), (46) equals the matrix product ${\Sigma }_{\mathit{vv}}{\Sigma }_{ll}^{-1}$, whereby the analogy to the GMM is stated again [34]. The proof that the trace of R equals the overall redundancy can be found in Appendix B.

3. Results

In this section the results of the trajectory shown in Section 2.2. will be analyzed by applying the simplified GHM (Table 1) on the approach described in Section 2.1. In a first step, the influence of the observation groups—accelerometer a, magnetometer m, estimated parameters of the previous epoch x, system noise c and gyroscope g—on the estimated orientation will be assessed by means of the r_i. Figure 3 shows a representative section of the calculated group redundancies for two different specifications of the sensor- and system noise standard deviations.

For the beginning the focus lies on the left part of Figure 3, where the original standard deviations are used (

Table 2. Standard deviations used for the sensor measurements and system noise.

Standard Deviations	Gyroscope	System Noise	Accelerometer	Magnetometer
Original	30°/s	10°/s	1 m/s2	5 μT
Adapted	60°/s	20°/s	0.5 m/s2	2.5 μT

). These are derived for each observation in the trajectory parts where the user walked straight. Using these standard deviations leads to a small partial redundancy of the observation groups (x, c, g) of the system equation in comparison to the ones of the measurement equation (a, m). This means that the estimated orientation mainly relies on the system equation. The reason for the high weight of the system equation is that the resulting trajectory should be smoothed [32]. The problem is that if the model assumptions made in the system equation do not capture the reality, the resulting deviations have high influence on the estimated orientation. There are mainly two possibilities to intervene in the partial redundancy respectively the inner reliability, which will be covered in the next two sections.

3.1. Adaption of the Stochastic Model

Adaption of the stochastic model means to change the standard deviations of the different observations. To align the group redundancies of the observation groups (Figure 3 left), the standard deviations of the sensors and system noise are changed until an improvement is visible. This leads to the group redundancies in the right part of Figure 3, where the adapted standard deviations of Table 2 were used. The weight of the system equation in comparison to the measurement equation is now reduced, whereas the r_i of the gyroscope and the system noise are now clearly higher. The r_i of the magnetometer get close to the ones of the gyroscope if it is used when the user is turning, which yields in good mutual controllability.

To identify the critical observations, the r_i of the individual observations are analyzed. Figure 4 shows, that the r_i of the accelerometer are continuously higher than 0.2, which means that controllability is sufficient. Calculating the thresholds from which deviations can be—statistically justified—detected (according to [29]), gives ~1.4 m/s² for the y-component of the accelerometer and ~2 m/s² for the x- and z-component. The user’s motion causes systematic deviations which exceeds these thresholds. Hence, accelerometer observations which should not be used for calculating pitch and roll can be detected.

The r_i of the magnetometer show a remarkable, alternating behavior. In the trajectory parts where the user walks straight, the r_i of m_y are close to zero. Whereas in parts where the user turns, controllability of m_x and m_z is bad. If there are systematic deviations in these observations, they cannot be detected and have a high influence on the estimated orientation angles. The r_i of x_ϕ and x_θ are again continuously higher than 0.2 and therefore sufficiently controlled. The same findings hold for c_ϕ and c_θ, whereas c_ψ is totally uncontrolled. Generally the r_i of the previously estimated parameters and the system noise behave in a similar manner. g_y is also uncontrolled as its r_i is very close to zero. The r_i of g_z are again higher than 0.2. The behavior of the r_i of the previously estimated yaw angle is interesting, as they go up to 0.25 if the gyroscope measurements are used in the system equation. If the gyroscope is not used in the following epochs the r_i decrease.

From the analysis of Figure 4 it can be concluded that there are uncontrollable observations in this approach, even if the standard deviations are adapted. To better understand the behavior of the r_i, the influence of the change of the standard deviation σ_i of one observation on the partial redundancy of the other observations is analyzed. This analysis should support an aimed change of the standard deviations to improve the r_i.

Table 3. Influence of the change of the standard deviation of one observation on the partial redundancy of the other observations. “+” indicates a raise and “−” a reduction of the partial redundancy. At the left sign, the standard deviation was increased by the factor 10 and at the right sign it was decreased by the factor 0.1. Grey shaded cells show the change in the corresponding observation and green means that the observations are functionally related.

	σg-y	σg-z	σc-ϕ	σc-θ	σc-ψ	σa-x	σa-y	σa-z	σm-x	σm-y	σm-z
rx-ϕ			−|−				−|+
rx-θ				−|−		−|+		−|+
rx-ψ		+|−			+|				−|+	−|+	−|+
rg-y	+|
rg-z		+|−							−|+	−|+	−|+
rc-ϕ			+|−				−|+			|+
rc-θ				+|−		−|+		−|+		|+
rc-ψ					+|−
ra-x				−|+		+|−		−|+
ra-y			−|+				+|−
ra-z				−|+		−|+		+|−		|+
rm-x		−|+							+|−	−|+	−|+
rm-y		−|+							−|+	+|−	−|+
rm-z		−|+							−|+	−|+	+|−

shows how the r_i change, if the standard deviation of one observation is multiplied with the factor 10 (left sign in Table 3) respectively 0.1 (right sign in Table 3). The standard deviation of all the observations was varied, except for the ones of the previously estimated parameters, as they are a direct result of the Kalman filter.

The grey shaded cells show the influence of a change of σ_i of one observation on the corresponding r_i. If σ_i is increased, the weight of the observation will be less in state estimation. Hence the corresponding r_i should also be increased. A reduction of σ_i should cause a smaller r_i on the contrary. This is the case for most of the observations, except for g_y. The reason therefore is that its r_i are already close to zero and a reduction of σ_i has no effect.

Green shaded cells show a direct functional relation of the corresponding observations (see (1)–(4) and (12)–(14)). The change of σ_i of one observation influences also the partial redundancy of other observations. A raise of one σ_i should cause a reduction of partial redundancy of other observations, as they get more weight in the state estimation. A reduction of one σ_i should raise the partial redundancy of other observations on the contrary. Though, the change of σ_i of c_ϕ and c_θ, show another behavior. By increasing as well as decreasing σ_c-ϕ and σ_c-θ, the r_x-ϕ and r_x-θ are reduced. A raise of σ_i of g_z respectively of c_ψ causes also a raise of the r_i of the previously estimated yaw angle.

A relation which is not expected from the system equations respectively the measurement equations, appears at the accelerometer measurements. a_y is only related to x_ϕ and c_ϕ, whereas a_x and a_z are related to x_θ and c_θ. The change of σ_i of one of the magnetometer measurements influences the r_i of x_ψ, g_z and the other two magnetometer components. m_y stands out, as a reduction of σ_i positively influences the r_i of c_ϕ and c_θ and the r_i of a_z. Hence, a reduction of σ_i of m_y would have the most positive influence on the r_i. As the controllability of this observation is bad, a further reduction is not advisable.

Changing the stochastic model clearly has an impact on the inner reliability but it also has limitations. As shown in Table 3, g_y as well as c_ψ are not controlled by any of the other observations (i.e., changing the standard deviation of any other observation does not influence r_g-y and r_c-ψ). The only way left in the stochastic model is to increase the standard deviation of such observations. In this case it is very likely, that the raised standard deviations cover the systematic deviations which actually should be detected.

3.2. Adaption of the Functional Model

In general, the functional model can be adapted by changing respectively extending the system equations and the measurement equations. As seen in the previous section, the y-component of the gyroscope and the system noise of the yaw angle stay problematical observations, as their controllability cannot be improved by changing σ_i of other observations. Through the example of the r_i of the gyroscope, the functional influences should be analyzed. (47) and (48) show the formulas to calculate the r_i of the gyroscope:

$$r_{g-y}={\sin }^2\phi \cdot \frac{{\sigma }_{g-y}^2\cdot d_{33}^{*}\cdot \Delta t^2}{{\cos }^2\theta }$$

(47)

$$r_{g-z}={\cos }^2\phi \cdot \frac{{\sigma }_{g-z}^2\cdot d_{33}^{*}\cdot \Delta t^2}{{\cos }^2\theta }$$

(48)

These equations are derived by evaluating (18), which is shown in Appendix C. d₃₃* corresponds to the third diagonal element of the inverse of D. The difference of (47) and (48) is the trigonometric function of the roll angle ϕ (σ²_g-y and σ²_g-z are assumed to be equal). As the roll angle is close to zero in the considered trajectory, also r_g-y will be close to zero even when changing σ_i of other observations (which influences d₃₃*). Hence, the smartphone orientation during the trajectory measurements has a huge impact on the calculated r_i—not only the ones from gyroscope, but also the ones from accelerometer and magnetometer ((12)–(14)). This seems reasonable if a closer look is taken on the quantities used to determine the smartphone orientation. The accelerometer should sense the gravity vector and the magnetometer the Earth magnetic field, which are both a vector quantity. The more sensor components sense these quantities, the better the mutual controllability should be. This should be the same for the gyroscope, which should sense the rotation of the user around the z-axis of the reference respectively navigation coordinate frame.

To check whether the assumptions made above are true, the recorded sensor data from the trajectory are rotated with the rotation matrix R given in (49):

$$\begin{array}{l}\mathit{R}={\mathit{R}}_x{\mathit{R}}_y{\mathit{R}}_z\\ {\mathit{R}}_x=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \phi & -\sin \phi \\ 0& \sin \phi & \cos \phi \end{array}\right], {\mathit{R}}_y=\left[\begin{array}{ccc}\cos \theta & 0& \sin \theta \\ 0& 1& 0\\ -\sin \theta & 0& \cos \theta \end{array}\right], {\mathit{R}}_z=\left[\begin{array}{ccc}\cos \psi & -\sin \psi & 0\\ \sin \psi & \cos \psi & 0\\ 0& 0& 1\end{array}\right]\end{array}$$

(49)

ϕ is chosen to be 45° and θ and ψ are 0°. Evaluating the rotated sensor data in the same Kalman filter as used above, gives the partial redundancies shown in Figure 5. The r_i of the previously estimated parameters, the system noise as well as the x-components of the sensors are not affected by the rotated data, when comparing the results with Figure 4. The y- and z-components of the sensors are more balanced now, leading to a better mutual controllability. Especially the mean level of the r_i of the gyroscope y-component are now approximately 0.1, which is clearly higher compared to Figure 4.

The idea of the following approach is to use multiple trajectory data, to determine the actual smartphone’s orientation. If one thinks of scenarios in crowded environments, it could be that multiple users have already taken the same path as the actual user does. If this data from the users who have taken the same path is stored, it could be used in a multiple trajectory data approach to determine the actual yaw angle, which equals the situation outlined in Section 2.3 where several sub-systems are used to determine the overall-system. Hence, the Gauss-Helmert Kalman filter which incorporates (32) in its functional model (33) has to be used now. In the case of orientation determination, the set of equations which unambiguously determine the state vector (31) are equal to (1), (2) and (3) respectively (4) for the actual trajectory. Now, for every additionally used trajectory, (1) and (2) have to be added to the unambiguous set of equations, as the smartphone’s orientation could be different (the only assumption is that there is the same path, e.g., the same yaw angle). Thus—if N_tr is the number of used trajectories—N_tr roll angles, N_tr pitch angles and one yaw angle have to be estimated. If there is a turn detected in the actual trajectory, gyroscope measurements from multiple smartphones have to be processed. Hence, every additional trajectory contributes with a formula of type (4) to the system equation, which leads to N_tr−1 condition Equations (32) in the system equation. The measurement equation consists of N_tr triples of (12)–(14).

In this article additional trajectories are simulated, which means that they are retrieved from the trajectory considered in this section by performing rotations in the same manner as done in (49). This is a simplification respectively a synthetical example whose aim is to retrieve more insights into how additional observations affect the partial redundancy. Before using additional rotated trajectories, the effect of using the same trajectory data two times will be analyzed. The reason is that no additional roll and pitch angles have to be estimated and therefore the results for the partial redundancy can be directly compared to the results from Section 3.1 respectively Figure 4. The results of the partial redundancy of the observations from the actual trajectory are shown in Figure 6.

The most obvious effect can be seen in the turn detection. The gyroscope is used less in the straight trajectory parts, which leads to slightly different appearance of the r_i. Nevertheless, the height of the r_i can be compared. The r_i of the previously estimated roll and pitch angel has not changed. It can be seen, that the r_i of the accelerometer, the system noise components and the previously estimated yaw angle are raised by using additional observations, which is especially important for the system noise of the yaw angle. The controllability of this quantity is now slightly improved. The components of the magnetometer and the gyroscope, which already had quite high r_i are now also raised, whereas an effect on the low r_i during turns is not visible.

This is not a surprise, as trajectory data is used two times where the controllability of these quantities is suboptimal. Figure 7 shows the partial redundancy of the observations from the actual trajectory by using one additional trajectory, where the smartphone’s orientation is different. This additional trajectory comes from rotating the trajectory data as mentioned above (49), which is a simplification respectively simulation. In a real-life application some sort of search-algorithm would have to be performed on the stored trajectories, to find the ones where users took the same path as in the actual trajectory. The only thing that has changed, are the r_i of the previously estimated roll and pitch angles, as they are clearly reduced now. The reason could be that the overall redundancy has not changed from the scenario in Figure 6 to the scenario in Figure 7 but more parameters have to be estimated. The remaining r_i have not changed. Especially r_g-y is still close to zero, despite using an additional trajectory where the smartphone’s orientation is different. The functional relation of the four gyroscope measurements—resulting from using the actual and the simulated trajectory (i.e., the two formulas of type (4))—doesn’t lead to an improved mutual controllability. The same holds for the r_i of the magnetometer measurements which are close to zero.

Figure 8 shows the ground truth from total station (see Section 2.2), the raw result from the magnetometer and different variants of the estimated yaw angle. For better visibility, again a representative section of the trajectory is chosen. Estimation variant 1 clearly differs from the rest of the estimation variants as it is much smoother. In this variant the original standard deviations from Table 2 are used, whereas in the other variants the adapted ones are used. The differences between variant 2 and 3 are negligible. This is reasonable, as the additional trajectory data of variant 3 is just rotated in comparison to the original trajectory data. Thus, the additional trajectory data used in this article does not contribute to the state estimation process but to the improvement of the inner reliability. In future experiments, real multiple trajectory data has to be collected from different smartphones and users to further evaluate the approach presented in this section. In this article, variant 1 performs slightly better with the drawback of worse inner reliability. In the beginning of the trajectory section shown in Figure 8, the raw result from magnetometer as well as the estimated yaw angles are shifted due to undetected systematic deviations. The aim of future research will be to detect such systematic deviations based on improved inner reliability, enabling better performance of a Kalman filter used for orientation determination.

4. Discussion

The main novelty of this article is the formulation of the Kalman filter and the redundancy matrix in the GHM. Depending on the system description in the system equation two different sets of update equations are derived, which are both applied to orientation determination as well as to partial redundancy calculation. The results of these derivations are used to analyze the inner reliability based on the partial redundancy. Analyzing the partial redundancy shows that in the Kalman filter used in Section 2.1 the system equation contribute disproportionately highly to the estimation of the orientation compared the measurement equation. This means that the observations of the system equation are nearly uncontrolled. Two general ways are considered to intervene into the inner reliability respectively partial redundancy.

First, effects of changes in the stochastic model are analyzed. By adapting the standard deviations, the group redundancies can be improved. Analyzing the partial redundancy of the individual observations shows that observations which crucially determine the orientation are still badly controlled. Systematic deviations in such observations can cause huge disturbances in the estimated quantities. By increasing and decreasing the standard deviations, their influence on the partial redundancy was analyzed. It appears that reducing the standard deviation of the magnetometer’s y-component, could increase many partial redundancies but its controllability is bad. The y-component of the gyroscope as well as the system noise of the yaw angle cannot be controlled by observations of the actual trajectory. It has to be mentioned that the resulting insights are in a way limited to the considered trajectory. Such adaptions in the stochastic model have to be also tested in other trajectories and different devices.

Furthermore, effects of changes in the functional model on the partial redundancy are considered. In a first step the influences of the smartphone’s orientation on the partial redundancy were analyzed. Rotating the trajectory data shows that if more than one sensor axis is sensing the quantity of interest, the controllability will be improved. Using multiple trajectory data should take advantage of this behavior. The formulation of the Kalman filter update equations in the GHM enables processing multiple trajectory data in a Kalman filter. Partial redundancy is improved through this approach but still there are critical observations, such as the y-components of the gyroscope and magnetometer. Just using more data does not mandatory lead to better controllability in the GHM. Additionally, improving inner reliability does not lead to better accuracy of the estimated quantities. These questions respectively challenges will be addressed in future studies.