Dual MIMU Pedestrian Navigation by Inequality Constraint Kalman Filtering

Abstract

The foot-mounted inertial navigation system is an important method of pedestrian navigation as it, in principle, does not rely any external assistance. A real-time range decomposition constraint method is proposed in this paper to combine the information of dual foot-mounted inertial navigation systems. It is well known that low-cost inertial pedestrian navigation aided with both ZUPT (zero velocity update) and the range decomposition constraint performs better than those in their own respective methods. This paper recommends that the separation distance between the position estimates of the two foot-mounted inertial navigation systems be restricted by an ellipsoidal constraint that relates to the maximum step length and the leg height. The performance of the proposed method is studied by utilizing experimental data, and the results indicate that the method can effectively correct the dual navigation systems’ position over the traditional spherical constraint.

1. Introduction

Positioning and tracking systems have developed over several decades in various applications, ranging from tracking of pedestrians to autonomous vehicles [1]. A high-precision navigation system is often needed for urban and indoor lives where GPS is unavailable, for example, in wearable body area networks [2]. Increasing attention for the pedestrian navigation problem can be partially attributed to the significant progress in affordable wearable computing platforms and enhancement in sensor quality, especially with respect to micro electro mechanical systems (MEMS) [3]. The micro inertial measurement units (MIMU) usually consist of gyroscopes, accelerometers, magnetometers, and pressure sensors [4].

In real-life application, however, the low-cost inertial navigation systems suffer from the accumulation of errors while calculating the traveled distance of the objects. These errors cause the trajectories to drift away from the actual path as time grows. An effective technique is to bind the error growth using ZUPT [5]. In [6], the authors presented an open-source, real-time, embedded implementation of a foot-mounted, zero-velocity-update-aided inertial navigation system (INS).

When we use the ZUPT-aided INS to track pedestrian location, the heading drift of navigation information is unobservable [7], so it is difficult to get accurate location information only through ZUPT in long-term navigation. Some previous works used the information fusion of multiple inertial sensors to ulteriorly correct the position [8,9]. As shown in [8], the use of the two-feet range constraints can significantly improve the navigation performance and a 110 m straight-line experiment showed that the spherical constraint algorithm can reduce the mean error and covariance of the final position estimates.

In [9], the authors proposed a method for a dual-mounted INS to reduce the systematic heading drift. They set up a dual MIMU system with two feet, in which the maximum distance between the two systems is derived from the calibration procedure. The proposed algorithm can obtain the walk trajectory when the initial heading estimates are known, but the method did not significantly improve the positioning accuracy of the system.

The experiments in [8,9] are two-dimensional (2-D) plane experiments, and the feasibility of the algorithm in three-dimensional space was not verified. When we use the spherical constraint [8] based on maximum step size to correct the position of the two feet, if the step size is greater than the height of heels, the spherical constraint algorithm would have little significant effect on the height constraint of the two feet. Therefore, in this paper we take into account different separation constraints in level and height directions, which is shown to be more effective to correct the pedestrian location.

2. Principle and Theory

2.1. Discrete Kalman Filter

The Kalman filter addresses the general problem of estimating the state $X_k$ of a discrete-time process that is governed by the linear stochastic difference equation [10,11]:

$$X_k={\Phi }_{k,k-1}{X}_{k-1}+{\Gamma }_{k-1}{W}_{k-1}$$

(1)

with a measurement $Z_k$ that is

$$Z_k=H_k{X}_k+V_k$$

(2)

where ${\Phi }_{k,k-1}$ denotes transition matrix relating the state at the previous time step $t_{k-1}$ to the state at the current step $t_k$, ${\Gamma }_{k-1}$ denotes the system noise drive matrix, $H_k$ represents a measurement matrix, $V_k$ is series of measurement noise, and $W_k$ represents noise excitation sequence for the system. The $W_k$ and $V_k$ simultaneously meet:

$$\begin{array}{l}E\left[W_k\right]=0,\mathrm{Cov}\left[W_k,W_j\right]=E\left[W_k{W}_j^T\right]=Q_k{\delta }_{kj}\\ E\left[V_k\right]=0,\mathrm{Cov}\left[V_k,V_j\right]=E\left[V_k{V}_j^T\right]=R_k{\delta }_{kj}\\ \mathrm{Cov}\left[W_k,V_j\right]=E\left[W_k{V}_j^T\right]=0\end{array}.$$

(3)

In practice, the process noise covariance matrices $Q_k$ and measurement noise covariance matrices $R_k$ might change with each time step or measurement, and we assume they are a positive definite matrix:

$${\widehat{X}}_k={\widehat{X}}_{k/k-1}+K_k\left(Z_k-H_k{\widehat{X}}_{k/k-1}\right)$$

(4)

$$K_k=P_{k/k-1}{H}_k^T{\left(H_k{P}_{k/k-1}{H}_k^T+R_k\right)}^{-1}$$

(5)

$$P_{k/k-1}={\Phi }_{k,k-1}{P}_{k-1}{\Phi }_{k.k-1}^T+{\Gamma }_{k-1}{Q}_{k-1}{\Gamma }_{k-1}^T$$

(6)

$$P_k=\left(I-K_k{H}_k\right)P_{k/k-1}$$

(7)

Equations (4)–(7) are the basic equations of Kalman filtering. If the initial values about ${\widehat{X}}_0$ and $P_0$ are given, we can perform the state estimation at time $k$ according to measurements $Z_k$ at the same moment.

2.2. Inequality Kalman Filter

The inequality Kalman filter appears in solving practical problems between state variables where there are inequality relationships [12]. The inequality relationship can be expressed as a constraint equation and combined with the Kalman filter, as a result of which the optimal solution strictly conforms to the inequality constraints between state variables, and a better result could be obtained.

The inequality-constrained Kalman optimal solution [13] is expressed as

$$\begin{array}{l}\underset{\widehat{x}}{\min }{\left({\widehat{x}}_k-x_k\right)}^{T}T\left({\widehat{x}}_k-x_k\right)\\ L{x}_k\le d\end{array}\}$$

(8)

where $x_k$ is the unconstrained (standard) Kalman filter estimate and $T$ is a symmetric positive definite weighting matrix [14,15]. As such, the weighted error of the constrained filter is minimized [16]. Expanding the first type of Equation (8):

$${\left({\widehat{x}}_k-x_k\right)}^{T}T\left({\widehat{x}}_k-x_k\right)={\widehat{x}}_k^{T}T{\widehat{x}}_k-2x_k^{T}T{\widehat{x}}_k+x_k^{T}T{x}_k$$

(9)

Thus, the inequality constrained problem can be further simplified as

$$\begin{array}{l}\underset{\widehat{x}}{\min }\left({\widehat{x}}_k^{T}T{\widehat{x}}_k-2x_k^{T}T{\widehat{x}}_k\right)\\ L{x}_k\le d\end{array}\}$$

(10)

3. Methods

3.1. Generalized Likelihood Ratio Test (GLRT)

The output of MIMU can be expressed as

$$x_k={\left[\begin{array}{cc}{x}_k^a& x_k^{\omega }\end{array}\right]}^T$$

where the specific force measurement vector $x_k^a\in {\mathsf{\Omega }}^3$ and the angular rate measurements vector $x_k^{\omega }\in {\mathsf{\Omega }}^3$. Assuming a series of measured value $y_n={\left\{x_k\right\}}_{\mathsf{k}=\mathsf{n}}^{\mathsf{n}+\mathsf{N}-1}$. We employ a double hypothesis testing as such, ${\mathsf{H}}_0:$ MIMU stationary, ${\mathsf{H}}_1:$ MIMU moving. The false alarm probability is expressed as

$${\mathsf{P}}_{\mathsf{F}\mathsf{A}}=\mathsf{P}\left\{{\mathsf{H}}_0|{\mathsf{H}}_1\right\}=\alpha$$

The detection probability is ${\mathsf{P}}_{\mathsf{D}}=\mathsf{P}\left\{{\mathsf{H}}_0|{\mathsf{H}}_0\right\}$. Two hypotheses' observation data probability density functions are, respectively, defined as $p\left(y_n;{\mathsf{H}}_0\right)$ and $p\left(y_n;{\mathsf{H}}_1\right)$.

The mathematical sensor model can be expressed as $x_k=s_k\left(\theta \right)+v_k$, where $s_k\left(\theta \right)={\left[\begin{array}{cc}{s}_k^a\left(\theta \right)& s_k^{\omega }\left(\theta \right)\end{array}\right]}^{\mathsf{T}}$ and $v_k={\left[\begin{array}{cc}{v}_k^a& v_k^{\omega }\end{array}\right]}^{\mathsf{T}}$, the force of MIMU is $s_k^a\left(\theta \right)\in {\mathsf{\Omega }}^3$, and MIMU angular rate is expressed as $s_k^{\omega }\left(\theta \right)\in {\mathsf{\Omega }}^3$. The symbol $\theta$ denotes the vector of unknown elements $v_k^a\in {\mathsf{\Omega }}^3$ accelerometers noise, $v_k^{\omega }\in {\mathsf{\Omega }}^3$ gyroscopes noise. Assume the noises follows zero mean Gaussian distribution, with noise covariance matrix $Z=E\left\{v_k{v}_k^{\mathsf{T}}\right\}=\left[\begin{array}{cc}{\sigma }_a^2{\mathsf{I}}_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& {\sigma }_{\omega }^2{\mathsf{I}}_{3\times 3}\end{array}\right]$, where ${\sigma }_a^2$ and ${\sigma }_{\omega }^2$, respectively, represent accelerometers and gyroscopes noise variance.

Since the sensor measurement can be obtained from the joint probability density as

$$p\left(y_n;\theta ,{\mathsf{H}}_i\right)=\prod _{k\in {\mathsf{\Omega }}_n}p\left(x_k^a;\theta ,{\mathsf{H}}_i\right)p\left(x_k^{\omega };\theta ,{\mathsf{H}}_i\right)$$

(11)

where:

$$p\left(x_k^a;\theta ,{\mathsf{H}}_i\right)=\frac{1}{{\left(2\pi {\sigma }_a^2\right)}^{3/2}}\exp \left\{-\frac{1}{2{\sigma }_a^2}{\Vert x_k^a-s_k^a\left(\theta \right)\Vert }^2\right\}$$

$$p\left(x_k^{\omega };\theta ,{\mathsf{H}}_i\right)=\frac{1}{{\left(2\pi {\sigma }_{\omega }^2\right)}^{3/2}}\exp \left\{-\frac{1}{2{\sigma }_{\omega }^2}{\Vert x_k^{\omega }-s_k^{\omega }\left(\theta \right)\Vert }^2\right\}$$

GLRT is determined by the hypothesis ${\mathsf{H}}_0$ if

$${\mathsf{L}}_G\left(y_n\right)=\frac{p\left(y_n;{\widehat{\theta }}_0,{\mathsf{H}}_0\right)}{p\left(y_n;{\widehat{\theta }}_1,{\mathsf{H}}_1\right)}>\lambda$$

(12)

where $\lambda$ denotes the threshold. In Equation (12), ${\widehat{\theta }}_0$ and ${\widehat{\theta }}_1$ represent the maximum likelihood estimate of the unknown element under the assumptions ${\mathsf{H}}_0$ and ${\mathsf{H}}_1$, respectively. Equation (12) can be simplified as

$$L_G\left(y_n\right)=\exp \left(-\frac{1}{2{\sigma }_a^2}\sum _{k\in {\mathsf{\Omega }}_n}{\Vert x_k^a-g\frac{{\overline{x}}_k^a}{\Vert {\overline{x}}_k^a\Vert }\Vert }^2-\frac{1}{2{\sigma }_{\omega }^2}\sum _{k\in {\mathsf{\Omega }}_n}{\Vert x_k^{\omega }\Vert }^2\right)$$

(13)

$$T\left(y_n\right)=-\frac{2}{N}\ln L_G\left(y_n\right)=\frac{1}{N}\sum _{k\in {\mathsf{\Omega }}_n}\left(\frac{1}{{\sigma }_a^2}{\Vert x_k^a-g\frac{{\overline{x}}_k^a}{\Vert {\overline{x}}_k^a\Vert }\Vert }^2+\frac{1}{{\sigma }_{\omega }^2}{\Vert x_k^{\omega }\Vert }^2\right).$$

(14)

$T\left(y_n\right)<\lambda$ means that the pedestrian is in a stationary state.

In practice, ZUPT can effectively aid inertial navigation system to remove long-time accumulated errors [5,17]. The velocity error of carrier is used as a concept [18,19]. When pedestrians stay static, the MIMU measured velocity is regarded as an error to correct the system using Kalman filtering.

The state error vector is defined as

$$\delta X={\left[\begin{array}{cc}\delta {\varphi }^{\mathsf{T}}& \begin{array}{cc}\begin{array}{cc}\delta {\omega }^{\mathsf{T}}& \delta r^{\mathsf{T}}\end{array}& \begin{array}{cc}\delta v^{\mathsf{T}}& \delta a^{\mathsf{T}}\end{array}\end{array}\end{array}\right]}^{\mathsf{T}}$$

which, respectively, represents the three-dimensional attitude error, gyro drift, position error, velocity error and accelerometer bias.

The zero-velocity correction Kalman filter model is

$$\{\begin{array}{l}\delta X_k={\Phi }_k\delta X_{k-1}+W_{k-1}\\ \delta Z_k=H_k\delta X_k+V_k\end{array}.$$

(15)

When the MIMU is stationary, the speed is zero, in theory; thus, the ZUPT speed measurement equation is

$$\delta Z_{v,k}=\mathsf{\Delta }{v}_k^b=v_k^b-{\left[\begin{array}{ccc}0& 0& 0\end{array}\right]}^{\mathsf{T}}$$

(16)

where the state transition matrix is given as

$${\Phi }_k=\left[\begin{array}{ccccc}{\mathsf{I}}_{3\times 3}& -\mathsf{\Delta }t{C}_{bk|k-1}^n& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& {\mathsf{I}}_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& {\mathsf{I}}_{3\times 3}& \mathsf{\Delta }t{\mathsf{I}}_{3\times 3}& 0_{3\times 3}\\ \mathsf{\Delta }tS\left(f_k^n\right)& 0_{3\times 3}& 0_{3\times 3}& {\mathsf{I}}_{3\times 3}& \mathsf{\Delta }t{C}_{bk|k-1}^n\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& {\mathsf{I}}_{3\times 3}\end{array}\right]$$

(17)

$$S\left(f_k^n\right)=\left[\begin{array}{ccc}0& -a_{zk}^n& a_{yk}^n\\ a_{zk}^n& 0& -a_{xk}^n\\ -a_{yk}^n& a_{xk}^n& 0\end{array}\right]$$

(18)

where $S\left(f_k^n\right)$ is the specific force anti-symmetric matrix, and $H_k=\left[\begin{array}{ccccc}{0}_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& {\mathsf{I}}_{3\times 3}& 0_{3\times 3}\end{array}\right]$.

3.2. The Ellipsoidal Constraint Method

Each foot are fixedly mounted by a MIMU. For regular human kinematics, the separation distance between the right and left feet cannot be larger than a quantity known as foot-to-foot maximum separation [8,9]. The maximum step size is a typical feature of pedestrian to walk and can be used to constrain the navigation error [20,21], namely, in addition to using ZUPT to improve the accuracy of pedestrian navigation. In specific, we decompose the constraint into three degrees of freedom and then use the obtained sub-constraints to correct the navigation system. Based on this intuition, we constrain the position estimate of right and left foot-mounted ZUPT-aided INSs.

According to the coordinate system identified of the MTI-G-700 units (3D motion tracking system, from Xsens Technologies B.V., Enschede, The Netherlands), the carrier coordinate system, as shown in Figure 1, shows the X_b axis is parallel to the surface of the MIMU, in the forward direction, and the Z_b axis is perpendicular to the MIMU surface, in the upward direction. In this dual-MIMU integrated navigation system, the navigation coordinate system’s X_n axis is forward, the Y_n axis points to the right, and the Z_n axis perpendicular to the X_nOY_n plane, upwards. The coordinates of the navigation subsystem bound to the feet are defined in the same way.

For two MIMU navigation systems, the $i=\mathsf{L},\mathsf{R}$, system real state is described as $x_k^i$ (including position, velocity, and attitude), the estimated state as ${\widehat{x}}_k^i$ at the time $k$, where $x_k^i\in {\mathsf{ℝ}}^{n_i},{\widehat{x}}_k^i\in {\mathsf{ℝ}}^{n_i}$.

The joint state vector is defined as

$$\{\begin{array}{l}{x}_k\stackrel{def}{=}{\left[\begin{array}{cc}{\left(x_k^{\mathsf{L}}\right)}^{\mathsf{T}}& {\left(x_k^{\mathsf{R}}\right)}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}\\ {\widehat{x}}_k\stackrel{def}{=}{\left[\begin{array}{cc}{\left({\widehat{x}}_k^{\mathsf{L}}\right)}^{\mathsf{T}}& {\left({\widehat{x}}_k^{\mathsf{R}}\right)}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}\end{array}$$

where ${\widehat{x}}_k\in {\mathsf{ℝ}}^m\left(n_1+n_2=m\right)$.

Letting the maximum step size of the pedestrian be given by $\gamma$, the real displacement difference between the two navigation systems should be less than or equal to $\gamma$. As the leg height is subject to certain constraints, during the pedestrian normal walking state, the positions of the right and left foot can be approximately constrained in a ellipsoid (Figure 2). The position of one foot is constrained within the circle of radius $\gamma$ in the XOY plane, and is confined within the circle with leg-related radius $h$ in XOZ and YOZ planes, both centered at the other foot (Figure 3).

Assuming that left foot is on the ground and the right foot is in movement at moment k (Figure 3), then we can calculate $\alpha$, defined as the angle between the position of the two feet in the XOY plane (in navigation coordinate system):

$$\alpha =\arctan \left|\frac{x_k^{\mathsf{R}}-x_k^{\mathsf{L}}}{y_k^{\mathsf{R}}-y_k^{\mathsf{L}}}\right|.$$

(19)

As we can see from Figure 2, there is a space azimuth $\beta$ between the right and the static left foot, we can calculate this angle by the positional relationship between the feet:

$$\beta =\arctan \frac{\left|z_k^{\mathsf{R}}-z_k^{\mathsf{L}}\right|}{\sqrt{{(y_k^{\mathsf{R}}-y_k^{\mathsf{L}})}^2+{(x_k^{\mathsf{R}}-x_k^{\mathsf{L}})}^2}}$$

(20)

Therefore, the ellipsoidal constraint correction algorithm between the feet can be defined as

$$\left[\begin{array}{c}{L}_s\\ L_h\end{array}\right]\cdot x_k\le \left(\begin{array}{c}\begin{array}{c}{\gamma }_x\\ {\gamma }_y\end{array}\\ h_z\end{array}\right)$$

(21)

Defining the matricies,

$$L_s=\left[\begin{array}{cccc}\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}& 0_{2\times 6}& \begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\end{array}& 0_{2\times 6}\end{array}\right]$$

$$L_h=\left[\begin{array}{cccccccccccccccccc}0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& -1& 0& 0& 0& 0& 0& 0\end{array}\right]$$

where ${\gamma }_x$ represents the real-time constraint value of the ellipsoid constraint on the X_n axis, so ${\gamma }_x=\gamma \sin \alpha$; ${\gamma }_y$ represents the real-time constraint value of the ellipsoid constraint on the Y_n axis, so ${\gamma }_y=\gamma \cos \alpha$; $h_z$ represents the real-time constraint value of the ellipsoid constraint on the Z_n axis, so $h_z=h\tan \beta$.

We assume that two navigation systems attitude is accurate in the current moment when the decompose step size constraint. When $\frac{{\Vert L_s\cdot x_k\Vert }^2}{{\gamma }^2}+\frac{{\Vert L_h\cdot x_k\Vert }^2}{h_z^2}>1$ can constraint the state to be satisfied with $\left\{x\in {ℝ}^m:\frac{{\Vert L_s\cdot x_k\Vert }^2}{{\gamma }^2}+\frac{{\Vert L_h\cdot x_k\Vert }^2}{{h_z}^2}\le 1\right\}$ the state modification is recommended as

$$\{\begin{array}{l}p\left({\widehat{x}}_k\right)\stackrel{def}{=}{\arg }_x\min {\left({\widehat{x}}_k-x\right)}^{\mathsf{T}}{P}_k^{-1}\left({\widehat{x}}_k-x\right)\\ \frac{{\Vert L_s\cdot x_k\Vert }^2}{{\gamma }^2}+\frac{{\Vert L_h\cdot x_k\Vert }^2}{h_z^2}\le 1\end{array}$$

(22)

where $P_k^{-1}$ denotes the Kalman filter estimated covariance matrix state.

Defining $L=\left[\begin{array}{cccc}\begin{array}{ccc}1/{\gamma }_x& 0& 0\\ 0& 1/{\gamma }_y& 0\\ 0& 0& 1/h_z\end{array}& 0_{3\times 6}& \begin{array}{ccc}-1/{\gamma }_x& 0& 0\\ 0& -1/{\gamma }_y& 0\\ 0& 0& -1/h_z\end{array}& 0_{3\times 6}\end{array}\right]$, Equation (22) can be written as

$$\{\begin{array}{l}p\left({\widehat{x}}_k\right)\stackrel{def}{=}{\arg }_x\min {\left({\widehat{x}}_k-x\right)}^{\mathsf{T}}{P}_k^{-1}\left({\widehat{x}}_k-x\right)\\ x_k^{\mathsf{T}}{L}^{\mathsf{T}}L{x}_k\le 1\end{array}.$$

(23)

The covariance matrix of the process measurement noise of the dual-MIMU integrated navigation system is

$$Q_k=\left[\begin{array}{cccc}{Q}_a& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& Q_{\omega }& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& Q_a& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& Q_{\omega }\end{array}\right]$$

(24)

where $Q_a=\left[\begin{array}{ccc}{\sigma }_{ax}^2& 0& 0\\ 0& {\sigma }_{ay}^2& 0\\ 0& 0& {\sigma }_{az}^2\end{array}\right]$ , $Q_{\omega }=\left[\begin{array}{ccc}{\sigma }_{\omega x}^2& 0& 0\\ 0& {\sigma }_{\omega y}^2& 0\\ 0& 0& {\sigma }_{\omega z}^2\end{array}\right]$ , ${\sigma }_a={\left[\begin{array}{cc}0.5& \begin{array}{cc}0.5& 0.5\end{array}\end{array}\right]}^{\mathsf{T}}$ , and ${\sigma }_{\omega }={\left[\begin{array}{cc}0.5& \begin{array}{cc}0.5& 0.5\end{array}\end{array}\right]}^{\mathsf{T}}\times \pi /180$.

The covariance matrix of the measurement noise of the dual-MIMU integrated navigation system is

$$R_k=\left[\begin{array}{cccccc}{\sigma }_{vx}^2& 0& 0& 0& 0& 0\\ 0& {\sigma }_{vy}^2& 0& 0& 0& 0\\ 0& 0& {\sigma }_{vz}^2& 0& 0& 0\\ 0& 0& 0& {\sigma }_{vx}^2& 0& 0\\ 0& 0& 0& 0& {\sigma }_{vy}^2& 0\\ 0& 0& 0& 0& 0& {\sigma }_{vz}^2\end{array}\right]$$

(25)

where ${\sigma }_v={\left[\begin{array}{cc}0.01& \begin{array}{cc}0.01& 0.01\end{array}\end{array}\right]}^{\mathsf{T}}$.

The sampling rate of the filter is 400 Hz.

4. Experiment

To compare the performance of the proposed algorithm with the existing in [8] (the spherical constraint method). Experiments are carried out using two MTI-G-700 units and the performance parameters of them are shown in

Table 1. The performance parameters of MTI-G-700.

Sensors	Accelerometer		Gyroscope
	Typ	Max	Typ	Max
Standard full range	50 $\mathsf{m}/{\mathsf{s}}^2$	-	450$°/\mathsf{s}$	-
Bias repeatability (1 year)	0.03 $\mathsf{m}/{\mathsf{s}}^2$	0.05 $\mathsf{m}/{\mathsf{s}}^2$	0.2$°/\mathsf{s}$	0.5$°/\mathsf{s}$
In-run bias stability	40 $\mathsf{\mu }\mathsf{g}$	-	10$°/\mathsf{h}$	-
Noise density	80 $\mathsf{\mu }\mathsf{g}/\sqrt{\mathsf{h}\mathsf{z}}$	150 $\mathsf{\mu }\mathsf{g}/\sqrt{\mathsf{h}\mathsf{z}}$	0.01$°/\mathsf{s}/\sqrt{\mathsf{h}\mathsf{z}}$	0.015$°/\mathsf{s}/\sqrt{\mathsf{h}\mathsf{z}}$
Non-linearity	0.03% FS	5% FS	0.01% FS	-

. The procedure is summarized as follows:

(1)

In a complex 2D environment: some closed trajectory containing a straight line path and turning eight times (turning angle: 90°).

(2)

In a complex 3D environment: a six-story staircase, and parts of corridors in the Sheng-Hua building at the Central South University. The walk strats at the first floor and ends at the sixth floor.

According to the experimenter gait characteristics, we set $\gamma =0.6$ m and $h=0.3$ m, respectively.

In the 2-D closed experiment, the ZUPT-aided INS can track the pedestrian feet positions (Figure 4), but the distance of two feet reaches about 6 m, which is unreasonable in pedestrian navigation. Comparing the three different sets of trajectories, the position estimation information under the ellipsoidal restriction can obtain the position estimates more accurately. Figure 5 shows the relative positions of the two feet in the $z$ axis direction. Since $h_z$ constrains the altitude difference, the feet height difference can be reduced near to the true value and the pedestrian location can be more accurate.

For estimation evaluation, we have chosen the root mean square error (RMSE) as an accuracy measure in this work. This is used to measure the difference between the actual values and the output of an estimator. For quantitative comparison, we only checked the starting and final positions, both in 2D and 3D, in all walking tests. We assume the starting position of the left foot and right foot as (0, 0.1, 0) and (0, –0.1, 0), hence, only the final estimates are inserted into the RMSE formula. Both of these quality indicators are given in comparison

Table 2. RMSE deduction by using the ellipsoidal constraint method for a closed path test.

Method	(L/R) 2D RMSE (m)	(L/R) 3D RMSE (m)	Remarks
Unconstraint	1.2640/0.9493	1.5473/1.2293	Time: 73 s Distance: 61.6 m Error rate(%): 0.93
Spherical constraint	0.6533/0.6194	0.8732/0.8482
Ellipsoidal constraint	0.5709/0.5953	0.6977/0.7174

. Note that the unconstrained method represented the results of only ZUPT corrections. Both spherical and ellipsoidal constraint methods can reduce the 2D and 3D trajectories error, but the latter obviously reduces the proportion more, and the correction effect is more obvious.

In the 3D upstairs experiment, the result shows the maximum step ellipsoidal constraint method reduces the error accumulation in the $z$-axis direction effectively. From the results presented in Figure 6 and Figure 7, we can observe that the proposed method in this paper can reduce the altitude difference of the feet position from 2.6 m to 0.56 m. Comparing to the spherical constraint method, the relative positions concentrate between –0.3 to 0.3 m which are obviously smaller than the spherical confinement results and are more suitable to the characteristics of the feet height on the stairs. By analyzing the results obtained in Figure 6 and Figure 7, we observe that the proposed algorithm can track the feet trajectory more accurately than the spherical constraint method.

In the indoor upstairs test, because of the objective factors of the irregular staircase, we are unable to accurately know the horizontal coordinate of the end point. However, the position in the Z_n axis can be accurately measured, so, in this experiment, we only analyze the root mean square error in the Z_n axis direction. The numerical results are given in

Table 3. RMSE deduction by using ellipsoidal constraint method for indoor upstairs test.

Method	Left-3D RMSE (m)	Right-3D RMSE (m)	Remarks
Unconstraint	1.6712	0.8249	Time: 157 s Height: 31.75 m Error rate(%): 1.71
Spherical constraint	0.7676	0.7874
Ellipsoidal constraint	0.6537	0.5414

.

5. Results

Low-cost inertial pedestrian navigation aided with both ZUPT and the range decomposition constraint performs better than those in their own respective method. In this paper, we decompose the maximum step length along the navigation coordinate axes in real time, and establish an ellipsoidal constraint more suitable for actual walking situations. Each sub-constraint changes along with different times and makes the aiding scheme of the step size more specific and accurate. Experimental tests on different paths show that the proposed ellipsoidal constraint method can effectively improve the position accuracy of pedestrian navigation.