A Novel Calibration Method for Gyro-Accelerometer Asynchronous Time in Foot-Mounted Pedestrian Navigation System

Abstract

Pedestrian Navigation System (PNS) is one of the research focuses of indoor positioning in GNSS-denied environments based on the MEMS Inertial Measurement Unit (MIMU). However, in the foot-mounted pedestrian navigation system with MIMU or mobile phone as the main carrier, it is difficult to make the sampling time of gyros and accelerometers completely synchronous. The gyro-accelerometer asynchronous time affects the positioning of PNS. To solve this problem, a new error model of gyro-accelerometer asynchronous time is built. The effect of gyro-accelerometer asynchronous time on pedestrian navigation is analyzed. A filtering model is designed to calibrate the gyro-accelerometer asynchronous time, and a zero-velocity detection method based on the rate of attitude change is proposed. The indoor experiment shows that the gyro-accelerometer asynchronous time is estimated effectively, and the positioning accuracy of PNS is improved by the proposed method after compensating for the errors caused by gyro-accelerometer asynchronous time.

1. Introduction

Outdoor pedestrian positioning mainly relies on the Global Navigation Satellite System (GNSS). In indoor GNSS-denied environments, a pedestrian navigation system (PNS) is an important positioning method to locate a user’s position [1,2]. There are two kinds of technical schemes for PNS: One is the active positioning method utilizing Ultra-Wideband (UWB) [3], Bluetooth [4], ZigBee [5], Wireless Fidelity (Wi-Fi) [6], etc. The other one is the passive positioning method utilizing the Strapdown Inertial Navigation System (SINS) based on three gyros and three accelerometers [7], which makes it realizable in terms of locating the user’s position without any other external devices. A common method for the passive positioning method is to fix the MEMS Inertial Measurement Unit (MIMU) on the foot, which is also called foot-mounted pedestrian navigation system [8,9].

Limited by the low precision of MEMS inertial devices, the velocity and position of PNS will diverge rapidly over time [10]. To avoid the accumulation of velocity errors and position errors, the zero velocity updates (ZUPT) algorithm has been proposed and widely used. The velocity is approximately equal to zero when the foot touches the ground, which can be employed in a Kalman filter to correct the errors of PNS. Zero-velocity detection is one of the main factors affecting the correction effect of the ZUPT algorithm [11]. In reference [12], Barsocchi analyzed the gait characteristics and succeeded in extracting temporal gait parameters based on a sensorized footwear which contains a 3-axis gyro, a 3-axis accelerometer, a 3-axis magnetometer and five pressure sensors. References [13,14,15] present four commonly used zero-velocity detection methods which are named the stance hypothesis optimal detector, the acceleration moving variance detector, the acceleration magnitude detector, and the angular rate energy detector. The models of the four zero-velocity detection methods are compared. The experimental results showed that the detection method which utilizes the synchronous signals of gyros and accelerometers has the highest performance, and the zero-velocity detection methods perform better at slow or normal gait speed.

In order to further improve the positioning accuracy of PNS, researchers have studied the application of other sensors and constraint algorithms in PNS. In reference [16], Ruiz et al. made use of the zero angular rate updates (ZARU) algorithm assisted by MIMU and RFID with extended Kalman filter, INS algorithm, and ZUPT algorithm as the framework. The experimental results showed that the positioning error is only 1.5 m in a long-term working. Xin et al. utilized magnetometers to correct the heading errors [17]. However, it’s unpractical since the magnetic field is complex and changeable in an indoor environment. Ji et al. proposed an adaptive heading correction algorithm for suppressing magnetic interference in SINS which can correct the heading error effectively in the environments with magnetic field interference [18]. Borenstein et al. proposed a method called heuristic drift reduction (HDR), which uses the building orientation to constrain the user’s trajectory [19]. Abdulrahim et al. proposed an algorithm for generating heading measurements from basic knowledge of the orientation of the building in which the pedestrian is walking, and the results showed that the absolute position error at the final position is below 5 m, about 0.1% of the total travelled distance in the 40-min experiment with a total distance of 3 km [20]. Tanigawa et al. utilized a barometer to achieve the high-fidelity height tracking. However, in fire and smoke environments, due to the thermal noise and quantization noise, the barometer is no longer credible and useable [21]. Gu et al. proposed an integrated method utilizing gait detection and hidden Markov height estimation algorithm, to realize the effective correction of height by taking advantage of the fixed height of floors and stairs without using barometers [22]. Jing et al. proposed an adaptive collaborative positioning algorithm which selects units for the collaborative network and integrates ranging measurement to constrain inertial measurement errors. The experimental results showed that the positioning accuracy is improved by 60% compared with the traditional method [23]. Ding and Skog et al. utilized MIMU arrays to improve the positioning accuracy [24,25]. However, these above methods do not consider the effect of the asynchronization of sampling time between gyros and accelerometers.

The traditional algorithms of SINS have been analyzed under the assumption that the sampling time between gyros and accelerometers is synchronous. However, in a practical system, the gyro signals need to pass a lowpass filter, and the frequency characteristics of gyros and accelerometers are inconsistent in the common case. Particularly, most of manufacturers which produce gyros or accelerometers seldom consider the compatibility between gyros and accelerometers. Therefore, the difference in phase-frequency characteristics will lead to the asynchronization of sampling time. In addition, most systems do not synchronize the sampling time of gyros and accelerometers by using the synchronous pulse. Therefore, the incompatibility between gyros and accelerometers and the circuit defect will cause the problem of gyro-accelerometer asynchronous time. If the sampling time of gyros and accelerometers is not synchronous, the velocity errors and position errors will be accumulated continuously over time. To solve this problem, Yan et al. analyzed the effect of gyro-accelerometer asynchronous time in SINS and successfully compensated for the velocity errors after calculating the parameter of gyro-accelerometer asynchronous time [26]. Wen et al. built a model of gyro-accelerometer asynchronous time in Dual-axis RINS, and proposed a calibration method that can improve the navigation velocity accuracy and system stability under long time navigation [27]. These experimental results showed that the methods proposed by Yan and Wen succeed in calibrating the parameter of gyro-accelerometer asynchronous time and compensating for the errors in a low-dynamic environment by making IMU rotate around a single axis or double axes regularly.

Nevertheless, the existing research work of gyro-accelerometer asynchronous time is not completely applicable in PNS. Compared with SINS, the rotation of MIMU in the foot-mounted pedestrian navigation system is more complex and intense. Such a high-dynamic environment will increase the effect of gyro-accelerometer asynchronous time on pedestrian navigation. Therefore, the traditional calibration methods and error models of gyro-accelerometer asynchronous time in SINS cannot be employed in PNS. However, there are few references of gyro-accelerometer asynchronous time in PNS. In the foot-mounted pedestrian navigation system with MIMU or mobile phone as the main carrier, the gyro-accelerometer asynchronous time is likely to reach tens of milliseconds, and the gyro-accelerometer asynchronous time is changeable each time the system is turned on. Due to the gyro-accelerometer asynchronous time, the velocity error tends to be nonlinear. Therefore, the traditional error models of PNS are no longer suitable. Besides, the traditional zero-velocity detection methods which utilize gyros and accelerometers do not consider the gyro-accelerometer asynchronous time. If the sampling time of gyros and accelerometers is asynchronous, the traditional zero-velocity detection methods may lead to missing detection or false detection.

In order to solve the above problems and improve the positioning accuracy of pedestrian navigation, the gyro-accelerometer asynchronous time in PNS is studied in this paper. Compared with the known and used methodologies, this paper focuses on the calibration method of gyro-accelerometer asynchronous time in PNS which is never studied before. The proposed technology intends to provide a new way to improve the positioning accuracy of foot-mounted pedestrian navigation system in which the communication infrastructure is damaged when fire, earthquake or other emergency occurs. The main contributions in this paper can be summarized as follows: (1) A new model of gyro-accelerometer asynchronous time is built, and the effect of gyro-accelerometer asynchronous time on pedestrian navigation is analyzed. (2) A calibration method of gyro-accelerometer asynchronous time and a new zero-velocity detection method based on the rate of attitude change are proposed.

The remainder of this paper is organized as follows: Section 2 introduces the reference frame definitions. Section 3 builds an error model of gyro-accelerometer asynchronous time and analyzes the effect of gyro-accelerometer asynchronous time on pedestrian navigation. Section 4 designs a Kalman filter to calibrate the gyro-accelerometer asynchronous time and proposes a zero-velocity detection method based on the rate of attitude change. Section 5 shows the experiments and results. The discussion and conclusion are presented in Section 6.

2. The Reference Frame Definitions

The reference coordinate frames used in this paper are defined as follows:

$o^i{x}^i{y}^i{z}^i$(i-frame): Earth-Centered Inertially Fixed (ECIF) orthogonal reference frame.

$o^e{x}^e{y}^e{z}^e$(e-frame): Earth-Centered Earth-Fixed (ECEF) orthogonal reference frame.

$o^b{x}^b{y}^b{z}^b$(b-frame): Body orthogonal reference frame aligned with Right-Forward-Up axes of MIMU.

$o^{b^{\prime }}{x}^{b^{\prime }}{y}^{b^{\prime }}{z}^{b^{\prime }}$($b^{\prime }$-frame): Accelerometer orthogonal reference frame aligned with accelerometer-sensitive axes.

$o^m{x}^m{y}^m{z}^m$(m-frame): Carrier orthogonal reference frame aligned with Right-Forward-Up axes of foot.

$o^h{x}^h{y}^h{z}^h$(h-frame): Horizontal orthogonal reference frame with X-axis and Y-axis in the local horizon.

$o^n{x}^n{y}^n{z}^n$(n-frame): Navigation orthogonal reference frame aligned with local East-North-Up geodetic axes.

3. Gyro-Accelerometer Asynchronous Time

3.1. Error Model of Gyro-Accelerometer Asynchronous Time

In the foot-mounted pedestrian navigation system, the MIMU rotates with the movement of foot. In theory, the gyro outputs and accelerometer outputs are both aligned with b-frame, which is calculated by gyros. However, due to the gyro-accelerometer asynchronous time, the accelerometer outputs are aligned with $b^{\prime }$-frame instead of b-frame. Therefore, the coordinate frames aligned with gyro-sensitive axes and accelerometer-sensitive axes are no longer consistent. The angle errors between b-frame and $b^{\prime }$-frame can be written as:

$$\left\{\begin{array}{c}\delta \theta =\theta -{\theta }^{\prime }\\ \delta \gamma =\gamma -{\gamma }^{\prime }\\ \delta \psi =\psi -{\psi }^{\prime }\end{array}\right.$$

(1)

where $\theta$, $\gamma$ and $\psi$ are the pitch angle, roll angle and yaw angle from b-frame to n-frame, ${\theta }^{\prime }$, ${\gamma }^{\prime }$ and ${\psi }^{\prime }$ are the pitch angle, roll angle and yaw angle from $b^{\prime }$-frame to n-frame.

Figure 1 shows the dynamic inconsistent error of the coordinate frame caused by gyro-accelerometer asynchronous time. If the gyro-accelerometer asynchronous time of three axes is different, $b^{\prime }$-frame is non-orthogonal. Since the three-axis MEMS accelerometer chip is used in PNS, it can be assumed that the three gyros and three accelerometers have the same characteristics of asynchronous time. That is, the three gyros and three accelerometers have the same asynchronous time.

In this paper, ${\mathit{f}}_{}^{b^{\prime }}$ denotes the accelerometer output vector in $b^{\prime }$-frame. ${\tilde{\mathit{f}}}_{}^b$ denotes the measured specific force vector. In view of gyro-accelerometer asynchronous time, the specific force measurements are the accelerometer outputs in $b^{\prime }$-frame. Therefore, the measured specific force vector can be expressed as:

$${\tilde{\mathit{f}}}_{}^b={\mathit{C}}_n^{b^{\prime }}{\mathit{f}}_{}^n$$

(2)

where ${\mathit{f}}_{}^n={\left[\begin{array}{ccc}{f}_E^n& f_N^n& f_U^n\end{array}\right]}^T$ is the projection of ${\mathit{f}}_{}^{b^{\prime }}$ onto n-frame, and ${\mathit{C}}_n^{b^{\prime }}$ is the transform matrix from n-frame to $b^{\prime }$-frame.

The measured specific force vector in n-frame can be written as:

$${\tilde{\mathit{f}}}_{}^n={\mathit{C}}_b^n{\tilde{\mathit{f}}}_{}^b$$

(3)

where ${\mathit{C}}_b^n$ is the transform matrix from b-frame to n-frame.

Therefore, the specific force error vector $\delta {\mathit{f}}_{}^n$ caused by gyro-accelerometer asynchronous time is:

$$\begin{array}{ccc}\hfill \delta {\mathit{f}}_{}^n& =& {\tilde{\mathit{f}}}_{}^n-{\mathit{f}}_{}^n\hfill \\ & \approx & {\mathit{C}}_b^n{\mathit{C}}_n^{b^{\prime }}{\mathit{f}}_{}^n-{\mathit{f}}_{}^n\hfill \\ & =& ({\mathit{C}}_b^n{\mathit{C}}_n^{b^{\prime }}-\mathit{I}){\mathit{f}}_{}^n\hfill \end{array}$$

(4)

where

$${\mathit{C}}_b^n=\left[\begin{array}{ccc}{\mathrm{c}}_{\psi }{\mathrm{c}}_{\gamma }-{\mathrm{s}}_{\psi }{s}_{\theta }{\mathrm{s}}_{\gamma }& -{\mathrm{s}}_{\psi }{\mathrm{c}}_{\theta }& {\mathrm{c}}_{\psi }{\mathrm{s}}_{\gamma }+{\mathrm{s}}_{\psi }{s}_{\theta }{\mathrm{c}}_{\gamma }\\ {\mathrm{s}}_{\psi }{\mathrm{c}}_{\gamma }+{\mathrm{c}}_{\psi }{\mathrm{s}}_{\theta }{\mathrm{s}}_{\gamma }& {\mathrm{c}}_{\psi }{\mathrm{c}}_{\theta }& {\mathrm{s}}_{\psi }{\mathrm{s}}_{\gamma }-{\mathrm{c}}_{\psi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\gamma }\\ -c_{\theta }{s}_{\gamma }& {\mathrm{s}}_{\theta }& {\mathrm{c}}_{\theta }{\mathrm{c}}_{\gamma }\end{array}\right],$$

$${\mathit{C}}_n^{b^{\prime }}=\left[\begin{array}{ccc}{\mathrm{c}}_{{\psi }^{\prime }}{\mathrm{c}}_{{\gamma }^{\prime }}-{\mathrm{s}}_{{\psi }^{\prime }}{s}_{{\theta }^{\prime }}{\mathrm{s}}_{{\gamma }^{\prime }}& {\mathrm{s}}_{{\psi }^{\prime }}{\mathrm{c}}_{{\gamma }^{\prime }}+{\mathrm{c}}_{{\psi }^{\prime }}{\mathrm{s}}_{{\theta }^{\prime }}{\mathrm{s}}_{{\gamma }^{\prime }}& -c_{{\theta }^{\prime }}{s}_{{\gamma }^{\prime }}\\ -{\mathrm{s}}_{{\psi }^{\prime }}{\mathrm{c}}_{{\theta }^{\prime }}& {\mathrm{c}}_{{\psi }^{\prime }}{\mathrm{c}}_{{\theta }^{\prime }}& {\mathrm{s}}_{{\theta }^{\prime }}\\ {\mathrm{c}}_{{\psi }^{\prime }}{\mathrm{s}}_{{\gamma }^{\prime }}+{\mathrm{s}}_{{\psi }^{\prime }}{s}_{{\theta }^{\prime }}{\mathrm{c}}_{{\gamma }^{\prime }}& {\mathrm{s}}_{{\psi }^{\prime }}{\mathrm{s}}_{{\gamma }^{\prime }}-{\mathrm{c}}_{{\psi }^{\prime }}{\mathrm{s}}_{{\theta }^{\prime }}{\mathrm{c}}_{{\gamma }^{\prime }}& {\mathrm{c}}_{{\theta }^{\prime }}{\mathrm{c}}_{{\gamma }^{\prime }}\end{array}\right],$$

$$\left\{\begin{array}{c}{\mathrm{s}}_{\alpha }=sin\left(\alpha \right)\hfill \\ {\mathrm{c}}_{\alpha }=cos\left(\alpha \right)\hfill \end{array}\right.(\alpha =\theta ,\gamma ,\psi ,{\theta }^{\prime },{\gamma }^{\prime },{\psi }^{\prime }).$$

In order to analyze the effect of gyro-accelerometer asynchronous time on pedestrian navigation conveniently and intuitively, the movement of foot is divided into three parts, including pitch motion, roll motion, and yaw motion. The effect of pitch motion on pedestrian navigation is analyzed as follows.

The gyro-accelerometer asynchronous time is denoted as $\tau$, and the initial three-axis attitude vector is denoted as ${\left[\begin{array}{c}0,0,{\psi }_0\end{array}\right]}^T$. If the MIMU rotates around the pitch axis, Equation (2) can be rewritten as:

$$\begin{array}{ccc}\hfill {\tilde{\mathit{f}}}_{}^b& =& {\mathit{C}}_n^{b^{\prime }}{\mathit{f}}_{}^n={\left({\mathit{C}}_{{\psi }_0}{\mathit{C}}_{{\theta }^{\prime }}\right)}^T{\mathit{f}}_{}^n\hfill \\ \hfill & =& \left[\begin{array}{ccc}1& 0& 0\\ 0& {\mathrm{c}}_{{\theta }^{\prime }}& {\mathrm{s}}_{{\theta }^{\prime }}\\ 0& -{\mathrm{s}}_{{\theta }^{\prime }}& {\mathrm{c}}_{{\theta }^{\prime }}\end{array}\right]\left[\begin{array}{ccc}{\mathrm{c}}_{{\psi }_0}& {\mathrm{s}}_{{\psi }_0}& 0\\ -{\mathrm{s}}_{{\psi }_0}& {\mathrm{c}}_{{\psi }_0}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{f}_E^n\\ f_N^n\\ f_U^n\end{array}\right]\hfill \\ \hfill & =& \left[\begin{array}{ccc}{\mathrm{c}}_{{\psi }_0}& {\mathrm{s}}_{{\psi }_0}& 0\\ -{\mathrm{s}}_{{\psi }_0}{\mathrm{c}}_{{\theta }^{\prime }}& {\mathrm{c}}_{{\psi }_0}{\mathrm{c}}_{{\theta }^{\prime }}& {\mathrm{s}}_{{\theta }^{\prime }}\\ {\mathrm{s}}_{{\psi }_0}{\mathrm{s}}_{{\theta }^{\prime }}& -{\mathrm{c}}_{{\psi }_0}{\mathrm{s}}_{{\theta }^{\prime }}& {\mathrm{c}}_{{\theta }^{\prime }}\end{array}\right]\left[\begin{array}{c}{f}_E^n\\ f_N^n\\ f_U^n\end{array}\right].\hfill \end{array}$$

(5)

Substituting Equation (5) into Equation (3), ${\tilde{\mathit{f}}}_{}^n$ can be expressed as:

$$\begin{array}{ccc}\hfill {\tilde{\mathit{f}}}_{}^n& =& {\mathit{C}}_b^n{\mathit{C}}_n^{b^{\prime }}{\mathit{f}}_{}^n={\mathit{C}}_{{\psi }_0}{\mathit{C}}_{\theta }{\mathit{C}}_n^{b^{\prime }}{\mathit{f}}_{}^n\hfill \\ \hfill & =& \left[\begin{array}{ccc}{\mathrm{c}}_{{\psi }_0}& -{\mathrm{s}}_{{\psi }_0}& 0\\ {\mathrm{s}}_{{\psi }_0}& {\mathrm{c}}_{{\psi }_0}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& {\mathrm{c}}_{\theta }& -{\mathrm{s}}_{\theta }\\ 0& {\mathrm{s}}_{\theta }& {\mathrm{c}}_{\theta }\end{array}\right]{\mathit{C}}_n^{b^{\prime }}{\mathit{f}}_{}^n\hfill \\ \hfill & \approx & \left[\begin{array}{ccc}1& 0& \delta \theta {\mathrm{s}}_{{\psi }_0}\\ 0& 1& -\delta \theta {\mathrm{c}}_{{\psi }_0}\\ -\delta \theta {\mathrm{s}}_{{\psi }_0}& \delta \theta {\mathrm{c}}_{{\psi }_0}& 1\end{array}\right]\left[\begin{array}{c}{\mathit{f}}_E^n\\ {\mathit{f}}_N^n\\ {\mathit{f}}_U^n\end{array}\right].\hfill \end{array}$$

(6)

For the convenience of analysis, it’s assumed that people walk, go upstairs or go downstairs at normal speed. Then $\delta \theta$, $\delta \gamma$ and $\delta \psi$ can be considered as small angles. Therefore, $\delta \theta$, $\delta \gamma$ and $\delta \psi$ can be written as:

$$\delta \alpha \approx {\omega }_{\alpha }\tau \begin{array}{cc}& (\end{array}\alpha =\theta ,\gamma ,\psi )$$

(7)

where ${\omega }_{\theta }$, ${\omega }_{\gamma }$ and ${\omega }_{\psi }$ are the angular velocities along pitch axis, roll axis, and azimuth axis respectively.

Substituting Equation (7) into Equation (6), ${\tilde{\mathit{f}}}_{}^n$ can be rewritten as:

$${\tilde{\mathit{f}}}_{}^n=\left[\begin{array}{ccc}1& 0& {\omega }_{\theta }\tau {\mathrm{s}}_{{\psi }_0}\\ 0& 1& -{\omega }_{\theta }\tau {\mathrm{c}}_{{\psi }_0}\\ -{\omega }_{\theta }\tau {\mathrm{s}}_{{\psi }_0}& {\omega }_{\theta }\tau {\mathrm{c}}_{{\psi }_0}& 1\end{array}\right]\left[\begin{array}{c}{f}_E^n\\ f_N^n\\ f_U^n\end{array}\right].$$

(8)

Compared with ${\mathit{f}}_{}^n$, the accelerometer measurement errors in n-frame can be expressed as:

$$\left[\begin{array}{c}\delta f_E^n\\ \delta f_N^n\\ \delta f_U^n\end{array}\right]=\left[\begin{array}{c}{\omega }_{\theta }\tau sin{\psi }_0{f}_U^n\\ -{\omega }_{\theta }\tau cos{\psi }_0{f}_U^n\\ -{\omega }_{\theta }\tau sin{\psi }_0{f}_E^n+{\omega }_{\theta }\tau cos{\psi }_0{f}_N^n\end{array}\right]$$

(9)

where $\delta {\mathit{f}}_{}^n={[\delta f_E^n,\delta f_N^n,\delta f_U^n]}^T$ denotes the accelerometer measurement error vector in n-frame.

According to [28], the velocity error model is:

$$\left\{\begin{array}{c}{\dot{\mathit{v}}}^n={\mathit{f}}_{}^n-(2{\omega }_{ie}^n+{\omega }_{en}^n)\times {\mathit{v}}^n+{\mathit{g}}^n\hfill \\ {\dot{\tilde{\mathit{v}}}}^n={\tilde{\mathit{f}}}_{}^n-(2{\omega }_{ie}^n+{\omega }_{en}^n)\times {\tilde{\mathit{v}}}^n+{\mathit{g}}^n\hfill \end{array}\right.$$

(10)

where ${\mathit{v}}^n={\left[\begin{array}{ccc}{v}_E^n& v_N^n& v_U^n\end{array}\right]}^T$ and ${\tilde{\mathit{v}}}^n={\left[\begin{array}{ccc}{\tilde{v}}_E^n& {\tilde{v}}_N^n& {\tilde{v}}_U^n\end{array}\right]}^T$ are the ideal velocity vector and calculated velocity vector in n-frame respectively, ${\mathit{g}}^n$ is the local gravity vector, ${\omega }_{ie}^n$ is the angular velocity vector of earth rotation and ${\omega }_{en}^n$ is the angular velocity vector of n-frame relative to earth.

Therefore, the velocity error vector caused by gyro-accelerometer asynchronous time can be expressed as:

$$\begin{array}{c}\hfill \delta {\mathit{v}}^n={\int }_0^t{\dot{\tilde{\mathit{v}}}}_{}^{n}dt-{\int }_0^t{\dot{\mathit{v}}}_{}^{n}dt\approx {\int }_0^t{\tilde{\mathit{f}}}_{}^n-{\mathit{f}}_{}^{n}dt={\int }_0^t\delta {\mathit{f}}_{}^{n}dt\end{array}$$

(11)

where $\delta {\mathit{v}}^n={\left[\begin{array}{ccc}\delta v_E^n& \delta v_N^n& \delta v_U^n\end{array}\right]}^T$ denotes the velocity error vector.

Substituting Equation (9) into Equation (11), the velocity error vector can be rewritten as:

$$\delta {\mathit{v}}^n=\tau {\int }_0^t{\omega }_{\theta }\left[\begin{array}{ccc}0& 0& sin{\psi }_0\\ 0& 0& -cos{\psi }_0\\ -sin{\psi }_0& cos{\psi }_0& 0\end{array}\right]\left[\begin{array}{c}{f}_E^n\\ f_N^n\\ f_U^n\end{array}\right]dt.$$

(12)

According to Equation (12), the horizontal velocity error is:

$$\begin{array}{c}\hfill \delta v_h^n=\sqrt{{\left(\delta v_E^n\right)}^2+{\left(\delta v_N^n\right)}^2}=\tau {\int }_0^t{\omega }_{\theta }{f}_U^{n}dt\end{array}$$

(13)

where the oritation is: $arctan(-\delta v_E^n/\phantom{\delta v_E^n\delta v_N^n}\delta v_N^n)={\psi }_0$.

The vertical velocity error is:

$$\begin{array}{c}\hfill \delta v_U^n=\tau {\int }_0^t{\omega }_{\theta }(-sin{\psi }_0{f}_E^n+cos{\psi }_0{f}_N^n)dt=\tau {\int }_0^t{\omega }_{\theta }{f}_{//}^{n}dt\end{array}$$

(14)

where $f_{//}^n$ is the equivalent velocity along the forward direction of the people’s movement. Similarly, $f_{\perp }^n$ denotes the equivalent velocity that is perpendicular to the forward direction of the people’s movement. $f_{\perp }^n$ can be expressed as:

$$f_{\perp }^n=cos{\psi }_0{f}_E^n+sin{\psi }_0{f}_N^n.$$

(15)

Further, the latitude error, longitude error and height error caused by gyro-accelerometer asynchronous time can be written as follows respectively:

$$\begin{array}{ccc}\hfill \delta L& =& {\int }_0^t\dot{\tilde{L}}dt-{\int }_0^t\dot{L}dt\hfill \\ & =& {\int }_0^t\frac{1}{R_M+\tilde{h}}{\tilde{v}}_N^{n}dt-{\int }_0^t\frac{1}{R_M+h}{v}_N^{n}dt\hfill \\ & \approx & -\frac{\tau cos{\psi }_0}{R_M+h}{\int }_0^t{\int }_0^t{\omega }_{\theta }{f}_U^{n}dtdt,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \delta \lambda & =& {\int }_0^t\dot{\tilde{\lambda }}dt-{\int }_0^t\dot{\lambda }dt\hfill \\ & =& {\int }_0^t\frac{secL}{R_N+\tilde{h}}{\tilde{v}}_E^{n}dt-{\int }_0^t\frac{secL}{R_N+h}{v}_E^{n}dt\hfill \\ & \approx & \frac{\tau sin{\psi }_{0}secL}{R_N+h}{\int }_0^t{\int }_0^t{\omega }_{\theta }{f}_U^{n}dtdt,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill \delta h& =& {\int }_0^t\dot{\tilde{h}}dt-{\int }_0^t\dot{h}dt\hfill \\ & =& {\int }_0^t{\tilde{v}}_U^{n}dt-{\int }_0^t{v}_U^{n}dt\hfill \\ & =& \tau {\int }_0^t{\int }_0^t{\omega }_{\theta }(-sin{\psi }_0{f}_E^n+cos{\psi }_0{f}_N^n)dtdt\hfill \end{array}$$

(18)

where L, $\lambda$, and h are the ideal latitude, longitude and height respectively, $\tilde{L}$, $\tilde{\lambda }$, and $\tilde{h}$ are the calculated latitude, longitude and height respectively, $R_M$ is the radius of curvature in meridian, and $R_N$ is the radius of curvature in prime vertical.

3.2. Effects of Gyro-Accelerometer Asynchronous Time on Pedestrian Navigation

According to Equation (12), the north velocity errors can be written as follows under the assumption of $\psi ={\psi }_0$ or $\psi ={\psi }_0+\pi$:

$$\left\{\begin{array}{c}{\left.\delta v_N^n\right|}_{\psi ={\psi }_0}=-\tau cos{\psi }_0{\int }_0^t{\omega }_{\theta }{f}_U^{n}dt\hfill \\ {\left.\delta v_N^n\right|}_{\psi ={\psi }_0+\pi }=\tau cos{\psi }_0{\int }_0^t{\omega }_{\theta }{f}_U^{n}dt.\hfill \end{array}\right.$$

(19)

Equation (19) shows that the north velocity errors are opposite in opposite directions. It’s indicated that the north velocity errors caused by gyro-accelerometer asynchronous time will cancel each other out in opposite directions. Further analysis shows that the gyro-accelerometer asynchronous time has the same effect on east velocity error.

According to Equation (13), the horizontal velocity error $\delta v_h^n$ is related to $\tau$, ${\omega }_{\theta }$ and $f_U^n$. Roughly, ${\omega }_{\theta }$ and $f_U^n$ can be regarded as sinusoidal signals approximately with the same period. Therefore, the horizontal velocity error will be accumulated, resulting in the drift of horizontal velocity in the forward direction of people’s movement.

Equation (14) shows that the vertical velocity error is related to $\tau$, ${\omega }_{\theta }$ and $f_{//}^n$. Since $f_{//}^n$ can be regarded as a sinusoidal signal approximately with the same period as ${\omega }_{\theta }$. Therefore, the vertical velocity error will be accumulated continuously, resulting in the rapid divergence of height eventually.

Similarly, when the MIMU rotates around the roll axis or azimuth axis, the effect of gyro-accelerometer asynchronous time on pedestrian navigation is analyzed and summarized respectively, and the results are shown in

Table 1. The effect of gyro-accelerometer asynchronous time on pedestrian navigation under different motions.

-	Pitch	Roll	Yaw
$\delta v_E^n$	$\tau sin{\psi }_0{\int }_0^t{\omega }_{\theta }{f}_U^{n}dt$	$\tau cos{\psi }_0{\int }_0^t{\omega }_{\gamma }{f}_U^{n}dt$	$-\tau {\int }_0^t{\omega }_{\psi }{f}_N^{n}dt$
$\delta v_N^n$	$-\tau cos{\psi }_0{\int }_0^t{\omega }_{\theta }{f}_U^{n}dt$	$\tau sin{\psi }_0{\int }_0^t{\omega }_{\gamma }{f}_U^{n}dt$	$\tau {\int }_0^t{\omega }_{\psi }{f}_E^{n}dt$
$\delta v_h^n$	$\tau {\int }_0^t{\omega }_{\theta }{f}_U^{n}dt$	$\tau {\int }_0^t{\omega }_{\gamma }{f}_U^{n}dt$	$\tau \sqrt{{\left({\int }_0^t{\omega }_{\psi }{f}_N^{n}dt\right)}^2+{\left({\int }_0^t{\omega }_{\psi }{f}_E^{n}dt\right)}^2}$
$\delta v_U^n$	$\tau {\int }_0^t{\omega }_{\theta }{f}_{//}^{n}dt$	$-\tau {\int }_0^t{\omega }_{\gamma }{f}_{\perp }^{n}dt$	0
$\delta L$	$-\frac{\tau cos{\psi }_0}{R_M+h}{\int }_0^t{\int }_0^t{\omega }_{\theta }{f}_U^{n}dtdt$	$\frac{\tau sin{\psi }_0}{R_M+h}{\int }_0^t{\int }_0^t{\omega }_{\gamma }{f}_U^{n}dtdt$	$\frac{\tau }{R_M+h}{\int }_0^t{\int }_0^t{\omega }_{\psi }{f}_E^{n}dtdt$
$\delta \lambda$	$\frac{\tau sin{\psi }_{0}secL}{R_N+h}{\int }_0^t{\int }_0^t{\omega }_{\theta }{f}_U^{n}dtdt$	$\frac{\tau cos{\psi }_{0}secL}{R_N+h}{\int }_0^t{\int }_0^t{\omega }_{\gamma }{f}_U^{n}dtdt$	$\frac{-\tau secL}{R_N+h}{\int }_0^t{\int }_0^t{\omega }_{\psi }{f}_N^{n}dtdt$
$\delta h$	$\tau {\int }_0^t{\int }_0^t{\omega }_{\theta }{f}_{//}^{n}dtdt$	$-\tau {\int }_0^t{\int }_0^t{\omega }_{\gamma }{f}_{\perp }^{n}dtdt$	0

.

As shown in Table 1, affected by gyro-accelerometer asynchronous time, when the MIMU rotates around the roll axis or azimuth axis, the north velocity errors or east velocity errors are opposite in opposite directions. It means that the horizontal position errors caused by gyro-accelerometer asynchronous time will cancel each other out in opposite directions. However, when the MIMU rotates around the roll axis, the vertical velocity error will be accumulated over time, leading to the rapid divergence of height eventually. The gyro-accelerometer asynchronous time makes no difference to the vertical velocity and height if the MIMU rotates around the azimuth axis.

In addition, Table 1 also shows the relationship between velocity errors and gyro-accelerometer asynchronous time, which provides a simple and convenient way to identify the parameter of gyro-accelerometer asynchronous time.

Table 2. The main characteristics of velocity errors under different motions.

Motion	Horizontal Velocity Error	Vertical Velocity Error
Pitch motion	The errors will cancel each other out in opposite directions	Increase
Roll motion	The errors will cancel each other out in opposite directions	Increase
Yaw motion	The errors will cancel each other out in opposite directions	Make no difference
Walking on flat ground	The errors will cancel each other out when walking on a closed-loop trajectory	Increase
Going upstairs	The errors will cancel each other out when walking on a closed-loop trajectory	Increase
Going downstairs	The errors will cancel each other out when walking on a closed-loop trajectory	Increase

shows the effect of gyro-accelerometer asynchronous time on velocity errors when the MIMU rotates around the pitch axis, roll axis, or azimuth axis. In practice, the MIMU rotates around three axes driven by foot. Therefore, it’s necessary to consider the effect of pitch motion, roll motion, and yaw motion on pedestrian navigation comprehensively under different gaits. Since the MIMU mainly rotates around the pitch axis, it’s reasonable to infer that the pitch motion is the main factor affecting velocity errors. Based on the above analysis, the major characteristics of velocity errors are summarized when people walk on flat ground, go upstairs, or go downstairs. The results are shown in Table 2. The horizontal velocity errors caused by gyro-accelerometer asynchronous time will cancel each other out when people walk on a closed-loop trajectory, while the horizontal velocity error will be accumulated with time when people walk on an open-loop trajectory, resulting in decreasing the horizontal positioning accuracy. The vertical velocity and height are mainly affected by pitch motion and roll motion. Since the MIMU rotates around the pitch axis and roll axis periodically, the vertical velocity error will be accumulated over time, leading to the rapid divergence of height eventually.

3.3. Simulation

In order to verify the correctness of the above analyses which are shown in Table 1 and Table 2, a simulated experiment was carried out. The gyro-accelerometer asynchronous time is set to 10ms, and the designed trajectory is divided into three sorts as follows:

1. People go up to the second floor from the start, then take a walk around the hall counterclockwise. 2. People continue to go up to the third floor from the second floor, then take a walk around the hall counterclockwise. 3. people continue to go up to the fourth floor from the third floor, then take a walk around the hall counterclockwise, and go down to the first floor from point A finally. The simulated trajectory is shown in Figure 2.

The yaw motion is modeled as a rotation around azimuth axis with an uniform angular velocity. The pitch motion and roll motion are modeled as periodic sinusoidal motions in one gait cycle:

$$\left\{\begin{array}{c}\theta =-{\theta }_{m}sin\left(2\pi f_{step}t\right)\hfill \\ \gamma ={\gamma }_{m}sin\left(2\pi f_{step}t\right)\hfill \end{array}\right.$$

(20)

where ${\theta }_m$ and ${\gamma }_m$ are the amplitudes of sinusoids, $f_{step}$ is the frequency corresponding to a gait cycle, and t is the working time.

Figure 3 presents the velocity errors and position errors caused by gyro-accelerometer asynchronous time, by comparing the ideal velocity (position) with the calculated velocity (position). It shows that the velocity errors have obvious characteristics in different positions when people walk, go upstairs, or go downstairs, which meet the main characteristics described in Table 2. In addition, Figure 3b shows that the gyro-accelerometer asynchronous time has a great influence on the position, which decreases the positioning accuracy of PNS. The rapid divergence of height makes it difficult to confirm exactly which floor the person is on. Therefore, it is of great significance to calibrate the gyro-accelerometer asynchronous time and compensate for the velocity errors and position errors for improving the positioning accuracy of PNS.

4. A Calibration Method for Gyro-Accelerometer Asynchronous Time

Considering that the gyro-accelerometer asynchronous time is changeable, a calibration method for gyro-accelerometer asynchronous time is proposed in this chapter.

4.1. Error Model of Pedestrian Navigation System Based on Gyro-Accelerometer Asynchronous Time

In theory, the projection of ${\mathit{f}}_{}^b$ onto n-frame can be written as ${\mathit{f}}_{}^n={\mathit{C}}_b^n{\mathit{f}}_{}^b$. Considering gyro-accelerometer asynchronous time, misalignment angles, and bias errors of accelerometers, the measured specific force vector in n-frame can be expressed as:

$${\tilde{\mathit{f}}}_{}^n={\tilde{\mathit{C}}}_b^n{\tilde{\mathit{f}}}_{}^{b^{\prime }}$$

(21)

where

$${\tilde{\mathit{f}}}_{}^{b^{\prime }}={\mathit{f}}_{}^{b^{\prime }}+\delta {\mathit{f}}_{}^b,{\tilde{\mathit{C}}}_b^n=[\mathit{I}-(\mathit{\varphi }\times )]{\mathit{C}}_b^n$$

where ${\delta \mathit{f}}^b$ is the bias error vector of accelerometers, $(·\times )$ denotes the operation of skew symmetric matrix, $\mathit{\varphi }$ is the attitude error vector.

According to [27], the attitude error model is:

$$\dot{\mathit{\varphi }}=\mathit{\varphi }\times \left({\omega }_{ie}^n+{\omega }_{en}^n\right)+{\delta \omega }_{ie}^n+{\delta \omega }_{en}^n-{\mathit{C}}_b^n{\delta \omega }_{ib}^b$$

(22)

where ${\delta \omega }_{ie}^n$,${\delta \omega }_{en}^n$ and ${\delta \omega }_{ib}^b$ are the error vectors of ${\omega }_{ie}^n$, ${\omega }_{en}^n$ and ${\omega }_{ib}^b$ respectively, and ${\omega }_{ib}^b$ is the gyro output vector.

Ignoring the second-order small errors, Equation (21) can be rewritten as:

$$\begin{array}{ccc}\hfill {\tilde{\mathit{f}}}_{}^n& =& [\mathit{I}-(\mathit{\varphi }\times )]{\mathit{C}}_b^n({\mathit{f}}_{}^{b^{\prime }}+\delta {\mathit{f}}_{}^b)\hfill \\ & =& [\mathit{I}-(\mathit{\varphi }\times )]{\mathit{C}}_b^n{\mathit{f}}_{}^{b^{\prime }}+[\mathit{I}-(\mathit{\varphi }\times )]{\mathit{C}}_b^n\delta {\mathit{f}}_{}^b\hfill \\ & \approx & [\mathit{I}-(\mathit{\varphi }\times )]{\mathit{C}}_{b^{\prime }}^n{\mathit{C}}_b^{b^{\prime }}{\mathit{f}}_{}^{b^{\prime }}+{\mathit{C}}_b^n\delta {\mathit{f}}_{}^b\hfill \end{array}$$

(23)

where ${\mathit{C}}_b^{b^{\prime }}$ is the transform matrix from b-frame to $b^{\prime }$-frame. The equivalent rotation vector from b-frame to $b^{\prime }$-frame is donated as $\delta \mathit{\vartheta }$ which can be expressed as:

$$\delta \mathit{\vartheta }\approx {\omega }_{nb}^b\tau \approx {\omega }_{ib}^b\tau$$

(24)

where ${\omega }_{nb}^b$ is the ideal angular velocity vector.

Therefore, ${\mathit{C}}_b^{b^{\prime }}$ can be formed as follows by using Rodrigues rotation formula [29]:

$$\begin{array}{ccc}\hfill {\mathit{C}}_b^{b^{\prime }}& =& {\mathit{M}}_{RV}\left(\delta \mathit{\vartheta }\right)\hfill \\ & =& \mathit{I}+\frac{sin\left(\left|\delta \mathit{\vartheta }\right|\right)}{\left|\delta \mathit{\vartheta }\right|}(\delta \mathit{\vartheta }\times )+\frac{1-cos\left(\left|\delta \mathit{\vartheta }\right|\right)}{{\left|\delta \mathit{\vartheta }\right|}^2}{(\delta \mathit{\vartheta }\times )}^2\hfill \\ & \approx & \mathit{I}+\frac{sin\left(\left|{\tilde{\omega }}_{nb}^b\tau \right|\right)}{\left|{\tilde{\omega }}_{nb}^b\tau \right|}[\left({\tilde{\omega }}_{nb}^b\tau \right)\times ]+\frac{1-cos\left(\left|{\tilde{\omega }}_{nb}^b\tau \right|\right)}{{\left|{\tilde{\omega }}_{nb}^b\tau \right|}^2}{[\left({\tilde{\omega }}_{nb}^b\tau \right)\times ]}^2\hfill \\ & \approx & \mathit{I}+\frac{sin\left(\left|{\omega }_{ib}^b\right|\tau \right)}{\left|{\omega }_{ib}^b\right|}({\omega }_{ib}^b\times )+\frac{1-cos\left(\left|{\omega }_{ib}^b\right|\tau \right)}{{\left|{\omega }_{ib}^b\right|}^2}{({\omega }_{ib}^b\times )}^2\hfill \end{array}$$

(25)

where ${\mathit{M}}_{RV}(·)$ is the function of equivalent rotation vector.

In a practical walking process, the angular velocity is several hundred degrees per second, while the gyro-accelerometer asynchronous time is several milliseconds or tens of milliseconds in most cases. So it’s reasonable to regard $\delta \mathit{\vartheta }$ as a small angle vector, and ${\mathit{C}}_b^{b^{\prime }}$ can be expressed as: ${\mathit{C}}_b^{b^{\prime }}\approx \mathit{I}+(\delta \mathit{\vartheta }\times )$. However, if people walk at high speed, such as running, sprinting, etc. The actual gyro outputs are likely to be more than one thousand degrees per second, then $\delta \mathit{\vartheta }$ is no longer a small angle vector, but a large angle vector. Therefore, it’s necessary to consider the nonlinear errors caused by gyro-accelerometer asynchronous time.

According to [30], the traditional velocity error model is:

$$\begin{array}{c}\hfill \delta {\dot{\mathit{v}}}^n=({\tilde{\mathit{f}}}_{}^n-{\mathit{f}}_{}^n)-(2{\tilde{\omega }}_{ie}^n+{\tilde{\omega }}_{en}^n)\times {\tilde{\mathit{v}}}^n+(2{\omega }_{ie}^n+{\omega }_{en}^n)\times {\mathit{v}}^n+({\tilde{\mathit{g}}}^n-{\mathit{g}}^n).\end{array}$$

(26)

Let:

$$\mathit{A}=\frac{sin\left(\left|{\omega }_{ib}^b\right|\tau \right)}{\left|{\omega }_{ib}^b\right|}({\omega }_{ib}^b\times ),$$

(27)

$$\mathit{B}=\frac{1-cos\left(\left|{\omega }_{ib}^b\right|\tau \right)}{{\left|{\omega }_{ib}^b\right|}^2}{({\omega }_{ib}^b\times )}^2.$$

(28)

Substituting Equations (23), (25), (27) and (28) into Equation (26), the velocity error model can be rewritten as:

$$\begin{array}{ccc}\hfill \delta {\dot{\mathit{v}}}^n& =& [\mathit{I}-(\mathit{\varphi }\times )]{\mathit{C}}_{b^{\prime }}^n{\mathit{C}}_b^{b^{\prime }}{\mathit{f}}_{}^{b^{\prime }}+{\mathit{C}}_b^n\delta {\mathit{f}}_{}^b-{\mathit{f}}_{}^n\hfill \\ & & -(2{\tilde{\omega }}_{ie}^n+{\tilde{\omega }}_{en}^n)\times {\tilde{\mathit{v}}}^n+(2{\omega }_{ie}^n+{\omega }_{en}^n)\times {\mathit{v}}^n+({\tilde{\mathit{g}}}^n-{\mathit{g}}^n)\hfill \\ & \approx & {\mathit{C}}_{b^{\prime }}^n(\mathit{I}+\mathit{A}+\mathit{B}){\mathit{f}}_{}^{b^{\prime }}-(\mathit{\varphi }\times ){\mathit{C}}_{b^{\prime }}^n(\mathit{I}+\mathit{A}+\mathit{B}){\mathit{f}}_{}^{b^{\prime }}+{\mathit{C}}_b^n\delta {\mathit{f}}_{}^b-{\mathit{f}}_{}^n\hfill \\ & & -(2{\tilde{\omega }}_{ie}^n+{\tilde{\omega }}_{en}^n)\times {\tilde{\mathit{v}}}^n+(2{\omega }_{ie}^n+{\omega }_{en}^n)\times {\mathit{v}}^n+({\tilde{\mathit{g}}}^n-{\mathit{g}}^n)\hfill \\ & \approx & \left[\left({\mathit{C}}_b^n{\mathit{f}}_{}^{b^{\prime }}\right)\times \right]\mathit{\varphi }+{\mathit{C}}_b^n(\mathit{A}+\mathit{B}){\mathit{f}}_{}^{b^{\prime }}+[{\mathit{C}}_b^n(\mathit{A}+\mathit{B}){\mathit{f}}_{}^{b^{\prime }}\times ]\mathit{\varphi }\hfill \\ & & -(2{\delta \omega }_{ie}^n+{\delta \omega }_{en}^n)\times {\mathit{v}}^n-(2{\omega }_{ie}^n+{\omega }_{en}^n)\times \delta {\mathit{v}}^n+{\mathit{C}}_b^n\delta {\mathit{f}}_{}^b+\delta {\mathit{g}}^n\hfill \end{array}$$

(29)

where:

$${\tilde{\omega }}_{ie}^n={\omega }_{ie}^n+{\delta \omega }_{ie}^n,{\tilde{\omega }}_{en}^n={\omega }_{en}^n+{\delta \omega }_{en}^n,$$

$${\tilde{\mathit{g}}}^n={\mathit{g}}^n+\delta {\mathit{g}}^n,{\tilde{\mathit{v}}}^n={\mathit{v}}^n+\delta {\mathit{v}}^n$$

where $\delta {\mathit{g}}^n$ is the error vector of ${\mathit{g}}^n$.

Equation (29) denotes the velocity error equation of PNS. The position errors can be modeled as:

$$\begin{array}{ccc}\hfill \delta \dot{L}& =& \dot{\tilde{L}}-L=\frac{1}{{\tilde{R}}_M+\tilde{h}}{\tilde{v}}_N^n-\frac{1}{R_M+h}{v}_N^n\hfill \\ & \approx & \frac{(v_N^n+\delta v_N^n)}{(R_M+\delta R_M)+(h+\delta h)}-\frac{1}{R_M+h}{v}_N^n,\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill \delta \dot{\lambda }& =& \dot{\tilde{\lambda }}-\dot{\lambda }=\frac{sec\tilde{L}}{{\tilde{R}}_N+\tilde{h}}{\tilde{v}}_E^n-\frac{secL}{R_N+h}{v}_E^n\hfill \\ & \approx & \frac{sec(L+\delta L)(v_E^n+\delta v_E^n)}{(R_N+\delta R_N)+(h+\delta h)}-\frac{secL}{R_N+h}{v}_E^n,\hfill \end{array}$$

(31)

$$\delta \dot{h}=\dot{\tilde{h}}-\dot{h}={\tilde{v}}_U^n-v_U^n=\delta v_U^n$$

(32)

where $\delta R_M$ and $\delta R_N$ are the errors of $R_M$ and $R_N$.

Since $R_M>>\delta L,\delta \lambda ,\delta h$, and $R_N>>\delta L,\delta \lambda ,\delta h$, Equations (30) and (31) can be rewritten as:

$$\delta \dot{L}\approx \frac{1}{R_M+h}\delta v_N^n-\frac{v_N^n}{{(R_M+h)}^2}\delta h,$$

(33)

$$\delta \dot{\lambda }\approx \frac{secL}{R_N+h}\delta v_E^n+\frac{v_E^{n}secLtanL}{R_N+h}\delta L-\frac{v_E^{n}secL}{{(R_N+h)}^2}\delta h.$$

(34)

4.2. Zero-Velocity Detection

Figure 4 presents the details of foot movement in a complete gait cycle, including four stages, named push-off phase, swing phase, hell strike phase, and stance phase [30]. The velocity is approximately equal to zero during the sole of foot touches the ground at the stance phase.

Reference [30] presents a zero-velocity detection method named Generalized Likelihood Ratio Test (GLRT) utilizing gyros and accelerometers comprehensively. However, GLRT is utilized under the assumption that there is no gyro-accelerometer asynchronous time. If people walk, especially at high speed, although the gyros measure the angular velocities at the stance phase, the accelerometers measure the specific force before it. Thus, it will lead to missing detection or false detection if gyro-accelerometer asynchronous time is ignored. GLRT performs better when people walk at slow or normal speed. However, if people walk at high speed, the performance of GLRT seems to be not enough. Considering the effect of gyro-accelerometer asynchronous time on zero-velocity detection, this paper proposes a zero-velocity detection method based on the rate of attitude change relying on gyros only, which is applicable under various gaits.

Compared with accelerometers, the gyro sampling signals hold the most reliable information for zero-velocity detection [12]. Therefore, the detection methods utilizing gyros only perform better than the detection methods utilizing accelerometers only. As shown in Figure 4, the pitch angle and roll angle are approximately equal to zero when the sole of foot touches the ground at the stance phase. Due to the low precision of MEMS inertial devices, the attitude errors will increase gradually with time. Therefore, it’s difficult to detect the zero-velocity based on attitudes. Since the change of pitch angle or roll angle is reliable in a short term, it’s an ideal way to detect the zero-velocity by utilizing the rates of change of pitch angle and roll angle.

Considering that there are installation angle errors when the MIMU is fixed on the foot, it is necessary to compensate for the installation errors in advance. Reference [31] presents a method to calibrate the pitch installation angle error and roll installation angle error by utilizing accelerometers. A self-calibration method of yaw installation angle error is shown as follows.

The installation angle error vector is denoted as ${\left[{\mathit{\varphi }}_x,{\mathit{\varphi }}_y,{\mathit{\varphi }}_z\right]}^T$. Using the chain rule of Direction Cosine Matrix (DCM) production, ${\mathit{C}}_h^n$ can be written as:

$${\mathit{C}}_h^n{=\mathit{C}}_m^n{\mathit{C}}_h^m.$$

(35)

Rewrite Equation (35) as:

$${\mathit{C}}_m^h={\left[{\left({\mathit{C}}_m^n\right)}^T{\mathit{C}}_h^n\right]}^T$$

(36)

where ${\mathit{C}}_m^n$ is the transform matrix from m-frame to n-frame, ${\mathit{C}}_h^n$ is the transform matrix from h-frame to n-frame, ${\mathit{C}}_m^h$ is the transform matrix from m-frame to h-frame. Using the chain rule of DCM production, ${\mathit{C}}_h^n$ can be expressed as:

$${\mathit{C}}_h^n={\mathit{C}}_b^n{\mathit{C}}_h^b={\mathit{C}}_b^n{\mathit{C}}_{{\mathit{\varphi }}_x}{\mathit{C}}_{{\mathit{\varphi }}_y}$$

(37)

where:

$${\mathit{C}}_{{\mathit{\varphi }}_x}=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\left({\mathit{\varphi }}_x\right)& sin\left({\mathit{\varphi }}_x\right)\\ 0& -sin\left({\mathit{\varphi }}_x\right)& cos\left({\mathit{\varphi }}_x\right)\end{array}\right],{\mathit{C}}_{{\mathit{\varphi }}_y}=\left[\begin{array}{ccc}cos\left({\mathit{\varphi }}_y\right)& 0& -sin\left({\mathit{\varphi }}_y\right)\\ 0& 1& 0\\ sin\left({\mathit{\varphi }}_y\right)& 0& cos\left({\mathit{\varphi }}_y\right)\end{array}\right].$$

${\mathit{C}}_m^h$ is related to ${\mathit{\varphi }}_z$ with the following relationship:

$$tan\left({\mathit{\varphi }}_z\right)=\frac{{\left({\mathit{C}}_m^h\right)}_{12}}{{\left({\mathit{C}}_m^h\right)}_{22}}.$$

(38)

Equations (36)–(38) present the calibration method of yaw installation angle error. However, it cannot be calculated if the foot is stationary. So it’s required to make MIMU rotate to realize the calculation of ${\mathit{\varphi }}_z$. This paper proposes a motion path of self-calibration for the yaw installation angle error, which is shown as follows.

It is assumed that the initial attitude vector is ${\left[\begin{array}{ccc}0& 0& {\phi }_0\end{array}\right]}^T$. When the person is walking along a straight line, the rotation angle vector from m-frame to n-frame is denoted as ${\mathit{\varphi }}_1\approx {\left[\begin{array}{ccc}{\theta }_1& {\gamma }_1& {\phi }_1\end{array}\right]}^T$ which represents the actual attitudes of the foot, and the rotation angle vector from h-frame to n-frame is denoted as ${\mathit{\varphi }}_2\approx {\left[\begin{array}{ccc}{\theta }_2& {\gamma }_2& {\phi }_2\end{array}\right]}^T$ which can be obtained from: ${\mathit{M}}_{RV}\left({\mathit{\varphi }}_2\right)={\mathit{C}}_b^n{\mathit{C}}_{{\mathit{\varphi }}_x}{\mathit{C}}_{{\mathit{\varphi }}_y}$. Since the person is walking along a straight line, the change of yaw angle can be regarded as a small angle. It’s reasonable to make an approximation that ${\theta }_1\approx {\theta }_2$, and ${\gamma }_1\approx {\gamma }_2$. Then, ${\mathit{\varphi }}_1$ can be expressed as:

$${\mathit{\varphi }}_1\approx {\left[\begin{array}{ccc}{\theta }_2& {\gamma }_2& {\phi }_0\end{array}\right]}^T.$$

(39)

Therefore, Equation (36) can be rewritten as:

$$\begin{array}{c}\hfill {\mathit{C}}_m^h\approx {\left\{{\left[{\mathit{M}}_{RV}\left({\mathit{\varphi }}_1\right)\right]}^T{\mathit{C}}_b^n{\mathit{C}}_{{\mathit{\varphi }}_x}{\mathit{C}}_{{\mathit{\varphi }}_y}\right\}}^T.\end{array}$$

(40)

Figure 5 presents the attitudes before and after compensating for the installation angle errors. It shows that the attitudes are corrected after compensating for the installation angle errors. The proposed calibration method of installation angle errors improves the recognition and characteristic of gait successfully.

Figure 6 presents the horizontal attitudes in a complete gait cycle. It’s indicated that the horizontal attitudes show different characteristics in different stages significantly. This paper proposes a zero-velocity detection method based on the rate of attitude change, which is shown as follows.

A fixed length sliding window is used, and the increments of pitch angle and roll angle of MIMU in $[t_k,t_{k+N}]$ are denoted as:

$$\left\{\begin{array}{c}\Delta {\theta }_i={\theta }_{i+1}-{\theta }_i\begin{array}{c}\hfill \end{array}(i=k,k+1,...,k+N)\hfill \\ \Delta {\gamma }_i={\gamma }_{i+1}-{\gamma }_i\begin{array}{c}\hfill \end{array}(i=k,k+1,...,k+N)\hfill \end{array}\right.$$

(41)

where N is the fixed length of sliding window.

If $\Delta {\theta }_i$ and $\Delta {\gamma }_i$ are less than the thresholds, it is indicated that the zero-velocity is detected. The judgment can be expressed as follows:

$$C_1=\left\{\begin{array}{c}1\begin{array}{c}\hfill \end{array}\forall t_i\in [t_k,t_{k+N}],\left|\Delta {\theta }_i\right|<G_{\theta }\hfill \\ 0\begin{array}{c}\hfill \end{array}others\hfill \end{array}\right.,$$

(42)

$$C_2=\left\{\begin{array}{c}1\begin{array}{c}\hfill \end{array}\forall t_i\in [t_k,t_{k+N}],\left|\Delta {\gamma }_i\right|<G_{\gamma }\hfill \\ 0\begin{array}{c}\hfill \end{array}others\hfill \end{array}\right.$$

(43)

where $G_{\theta }$ and $G_{\gamma }$ are the thresholds.

If $C_1\&C_2=1$, the zero-velocity is detected.

Figure 7 presents the results of zero-velocity detection under different gaits. The zero value of the purple line means that the foot is stationary at that moment, while the non-zero value of the purple line means that the foot is non-stationary at that moment. It shows that the zero-velocity is detected effectively by utilizing the proposed zero-velocity detection method.

4.3. Kalman Filter Design

Equation (29) shows that the velocity error model is nonlinear. Therefore, an Extended Kalman Filter (EKF) is utilized to estimate the parameter of gyro-accelerometer asynchronous time in this paper. The integrated framework block diagram of PNS is shown in Figure 8.

The state vector for Kalman filter is:

$$\mathit{X}={\left[\begin{array}{cccccc}{\mathit{\varphi }}^{\mathrm{T}}& {\left(\delta {\mathit{v}}^n\right)}^{\mathrm{T}}& {\left({\delta \mathit{p}}^n\right)}^{\mathrm{T}}& {\left({\epsilon }^b\right)}^{\mathrm{T}}& {\left({\mathbf{\nabla }}^b\right)}^{\mathrm{T}}& \tau \end{array}\right]}^{\mathrm{T}}$$

(44)

where ${\left(\delta {\mathit{p}}^n\right)}^{\mathrm{T}}={\left[\begin{array}{ccc}\delta L& \delta \lambda & \delta h\end{array}\right]}^T$ is the position error vector, ${\left({\epsilon }^b\right)}^{\mathrm{T}}={\left[\begin{array}{ccc}{\epsilon }_x^b& {\epsilon }_y^b& {\epsilon }_z^b\end{array}\right]}^T$ is the bias error vector of gyros, and ${\left({\mathbf{\nabla }}^b\right)}^{\mathrm{T}}={\left[\begin{array}{ccc}{\nabla }_x^b& {\nabla }_y^b& {\nabla }_z^b\end{array}\right]}^T$ is the bias error vector of accelerometers. ${\left({\epsilon }^b\right)}^{\mathrm{T}}$, ${\left({\mathbf{\nabla }}^b\right)}^{\mathrm{T}}$ and $\tau$ have the following relationship:

$$\left\{\begin{array}{c}{\dot{\epsilon }}^b=0\hfill \\ {\dot{\mathbf{\nabla }}}^b=0\hfill \\ \dot{\tau }=0\hfill \end{array}\right..$$

(45)

The state transition equation can be expressed as:

$$\dot{\mathit{X}}=\mathit{f}\left(\mathit{X}\right)+\mathit{GW}$$

(46)

where $\mathit{f}\left(\mathit{X}\right)$ is nonlinear vector function, $\mathit{G}$ is process noise coupling matrixr, $\mathit{W}$ is process noise vector.

The measurement equation can be written as:

$$\mathit{Z}=\mathit{HX}+\mathit{V}$$

(47)

where $\mathit{V}$ is measurement noise vector, $\mathit{Z}$ and $\mathit{H}$ are measurement vector and measurement matrix respectively.

In this paper, the height estimation method proposed by Gu et al. is utilized to estimate accurate height measurements for EKF [22]. Therefore, $\mathit{Z}$ and $\mathit{H}$ can be expressed as:

$$\mathit{Z}=\left[\begin{array}{c}{\tilde{\mathit{v}}}_{}^n-0\hfill \\ {\tilde{h}}_{}-h_{building}\hfill \end{array}\right],$$

(48)

$$\mathit{H}=\left[\begin{array}{ccccc}{\mathbf{0}}_{3\times 3}\hfill & {\mathit{I}}_{3\times 3}& {\mathbf{0}}_{3\times 2}& {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 7}\\ {\mathbf{0}}_{1\times 3}\hfill & {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{3\times 2}& 1& {\mathbf{0}}_{1\times 7}\end{array}\right]$$

(49)

where $h_{building}$ is the height measurement calculated by the height estimation method.

Discretizing Equations (46) and (47), the state transition equation and measurement equation can be rewritten as:

$$\left\{\begin{array}{c}{\mathit{X}}_k=\mathit{f}\left({\mathit{X}}_{k-1}\right)+{\Gamma }_{k-1}{\mathit{W}}_{k-1}\hfill \\ {\mathit{Z}}_k={\mathit{H}}_k{\mathit{X}}_k+{\mathit{V}}_k\hfill \end{array}\right.$$

(50)

where ${\mathit{W}}_{k-1}$ and ${\mathit{V}}_k$ are zero-mean Gaussian white noise vector sequences, both of which are not correlated. The feedback compensation of EKF can be expressed as:

$$\left\{\begin{array}{c}{\mathit{C}}_b^n=[{\mathit{I}}_{3\times 3}+(\mathit{\varphi }\times )]{\tilde{\mathit{C}}}_b^n\hfill \\ {\mathit{v}}^n={\tilde{\mathit{v}}}^n-\delta {\mathit{v}}^n\hfill \\ h=\tilde{h}-\delta h\hfill \end{array}\right..$$

(51)

5. Experiments and Analysis

The experiment was conducted in a building in Jimei District, Xiamen. As shown in Figure 9, to collect the data, the MIMU was fixed on the instep of the right foot. Then the experimenter walked on a designed trajectory that contains different gaits, including walking on flat ground, going upstairs, and going downstairs. While the experimenter was walking indoor, the MIMU sent the gyro outputs and accelerometer outputs to the Bluetooth receiver connected to a PC through the Bluetooth transmitter. In this experiment, a self-designed MIMU which contains three gyros and three accelerometers with large ranges is used. The performance of the MIMU is shown in

Table 3. Performance of MIMU.

Performance	Gyros	Accelerometers
In-run stability	${10}^{\circ }$/h	40 ug
Random walk	$0.4^{\circ }$/$\sqrt{}h$	0.06 m/s /$\sqrt{}h$
Full range	$\pm {2000}^{\circ }$/h	±40 g

.

The experimental environment is shown in Figure 10. The height of the floors is 4 m, and a start point and three standard points are set up in advance according to the engineering drawing of the building. Their relative coordinates in the building can be expressed as follows:

$P_0(0,0,0)$: the start point on the first floor; $P_1(0,0,8)$: standard point 1 on the third floor; $P_2(0,0,12)$: standard point 2 on the fourth floor; $P_3(0,0,8)$: standard point 3 on the third floor coincided with standard point 1.

The designed trajectory can be summarized as follows: The experimenter went up to the third floor from $P_0$, and took a walk around the hall counterclockwise to reach $P_1$. Then the experimenter went up to the fourth floor, and took a walk around the hall clockwise to reach $P_2$. Finally, the experimenter went down to the third floor, and took a walk around the hall counterclockwise to reach $P_3$. The total travelled distance is about 300 m.

After data collection is completed, the gyro sampling signals and accelerometer sampling signals are compared offline. The estimation of gyro-accelerometer asynchronous time is about 10 ms. A method named interpolation is utilized to compensate for the errors caused by gyro-accelerometer asynchronous time then. The error curves are drawn in Figure 11, by comparing the compensated velocity (position) with the uncompensated velocity (position). Figure 11 shows that the characteristics of velocity errors and position errors caused by gyro-accelerometer asynchronous time meet the results described in Table 1 and Table 2. It further proves the correctness of the above analyses.

Reference [32] presents the ZUPT algorithm and the height constraint algorithm, both of which show excellent performance in PNS. To evaluate the performance of the proposed calibration method, four methods are designed in this paper. Method 1 utilizes the ZUPT algorithm with gyro-accelerometer asynchronous ignored. Method 3 utilizes the ZUPT algorithm and the height constraint algorithm with gyro-accelerometer asynchronous ignored. Method 2 and method 4 are designed to verify the effectiveness and feasibility of the proposed calibration method by comparing the positioning accuracy with method 1 and method 3. The designed methods are shown in

Table 4. Experimental methods.

Method	Detail
Method 1	ZUPT with gyro-accelerometer asynchronous time ignored
Method 2	ZUPT with gyro-accelerometer asynchronous time considered
Method 3	ZUPT/height constraint with gyro-accelerometer asynchronous time ignored
Method 4	ZUPT/height constraint with gyro-accelerometer asynchronous time considered

in detail.

Figure 12 presents the estimations of gyro-accelerometer asynchronous time with method 2 and method 4. Since method 1 and method 3 do not consider the gyro-accelerometer asynchronous time, the information about method 1 or method 3 is not presented in Figure 12. It shows that the estimation of gyro-accelerometer asynchronous time with method 2 is about 8.92 ms, and the estimation of gyro-accelerometer asynchronous time with method 4 is about 8.31 ms. These results are very close to the offline calculated estimation. Moreover, Figure 12 also shows that the curves are convergent after about 60 s due to the complex movement of MIMU driven by foot, indicating that the gyro-accelerometer asynchronous time can be inspired in a relatively short time.

Figure 13 presents the calculated trajectories with method 1 and method 2, and Figure 14 presents the calculated trajectories with method 3 and method 4. Figure 13a shows that the height error with method 1 is positive while the height error with method 2 is negative. Further analysis shows that the height error with method 2 is affected by the low precision of MEMS inertial devices and the computational errors of inertial navigation algorithm, while the height error with method 1 is also affected by the gyro-accelerometer asynchronous time. Therefore, the height errors with method 1 and method 2 are opposite. Figure 13b shows that the calculated trajectory with method 2 is more consistent with the real trajectory than the calculated trajectory with method 1. In addition, as shown in Figure 14, compared with the calculated trajectories with method 1 and method 2, method 3 and method 4 have a higher positioning accuracy. The height constraint algorithm makes it possible to confirm exactly which floor the experimenter is on.

In order to further verify the effectiveness of the proposed calibration method, the positioning errors are compared. Figure 15 presents the position curves by four methods respectively. Figure 16 presents the position errors by four methods at three standard points. The details of the position errors are shown in

Table 5. Comparison of position errors by different methods ($P_1$).

Errors	Method 1	Method 2	Method 3	Method 4
North position error (m)	0.38	0.35	0.17	0.18
East position error (m)	0.36	0.23	0.27	0.15
Horizontal position error (m)	0.52	0.42	0.32	0.23
Error percentage (%D)	0.50	0.40	0.31	0.22
Height error (m)	1.03	0.76	0.08	0.06

,

Table 6. Comparison of position errors by different methods ($P_2$).

Errors	Method 1	Method 2	Method 3	Method 4
North position error (m)	0.93	0.75	0.44	0.36
East position error (m)	1.06	0.92	1.01	0.85
Horizontal position error (m)	1.41	1.19	1.10	0.92
Error percentage (%D)	0.70	0.59	0.54	0.46
Height error (m)	2.48	1.36	0.20	0.18

and

Table 7. Comparison of position errors by different methods ($P_3$).

Errors	Method 1	Method 2	Method 3	Method 4
North position error (m)	1.70	1.38	0.89	0.66
East position error (m)	1.84	1.60	1.66	1.43
Horizontal position error (m)	2.51	2.11	1.88	1.58
Error percentage (%D)	0.84	0.70	0.64	0.53
Height error (m)	3.37	2.11	0.22	0.23

. The experimental results show that the horizontal position errors with method 1 are 0.52 m, 1.41 m and 2.51 m, accounting for 0.50%, 0.70% and 0.84% of the total travelled distance respectively, and the horizontal position errors with method 2 are 0.42 m, 1.19 m and 2.11 m, accounting for 0.40%, 0.59% and 0.70% of the total travelled distance respectively. Compared with method 1, the horizontal positioning accuracy with method 2 increases by 23.81%, 18.49% and 18.96% respectively. The horizontal position errors with method 3 are 0.32 m, 1.10 m and 1.88 m, accounting for 0.31%, 0.54% and 0.64% of the total travelled distance respectively, and the horizontal position errors with method 4 are 0.23 m, 0.92 m and 1.58 m, accounting for 0.22%, 0.46% and 0.53% of the total travelled distance respectively. Compared with method 3, The horizontal positioning accuracy with method 4 increases by 39.13%, 19.57% and 18.99% respectively. Therefore, the results prove that the horizontal positioning accuracy is improved after compensating for the errors caused by gyro-accelerometer asynchronous time. Further analysis shows that the velocity errors are only corrected at the stance phase by using the ZUPT algorithm. Since the zero-velocity interval accounts for a small part of a complete gait cycle, the velocity errors and position errors will be accumulated continuously. Although the horizontal velocity errors caused by gyro-accelerometer asynchronous time will cancel each other out partly on a closed-loop trajectory, the complex movement of foot will still lead to the drift of velocity. Therefore, the horizontal positioning accuracy can be improved by compensating for the errors caused by gyro-accelerometer asynchronous time.

Besides, compared with method 1, the vertical positioning accuracy with method 2 increases by 0.27 m, 1.12 m and 1.26 m respectively at $P_0$, $P_1$ and $P_2$, indicating that the proposed calibration method can effectively compensate for the height error caused by gyro-accelerometer asynchronous time when using the ZUPT algorithm. However, the improvement of vertical positioning accuracy with method 4 is limited, due to the effective correction of height when using the height constraint algorithm.

6. Discussion and Conclusions

The experimental results show that the gyro-accelerometer asynchronous time can be estimated in a short time due to the high-dynamic environments. The proposed calibration method performs well in estimating the parameter of gyro-accelerometer asynchronous time. The positioning accuracy can be improved effectively after compensating for the errors caused by gyro-accelerometer asynchronous time. Since the position errors caused by gyro-accelerometer asynchronous time will be accumulated over time, the proposed calibration method contributes to improving the stability of PNS in a long time.

In the foot-mounted pedestrian navigation system with MIMU or mobile phone as the main carrier, the difference of phase-frequency characteristics between gyros and accelerometers will lead to the asynchronization of sampling time, resulting in the accumulation of velocity errors and position errors. To solve this problem, in this paper, an error model of gyro-accelerometer asynchronous time is built. The effect of gyro-accelerometer asynchronous time on pedestrian navigation and the main characteristics of velocity errors under different motions are analyzed. A filtering model is designed to calibrate the gyro-accelerometer asynchronous time via a Kalman filter. To avoid the missing detection and false detection, a zero-velocity detection method based on the rate of attitude change is proposed. The results of the 300 m-long experiment show that the gyro-accelerometer asynchronous time is estimated effectively and the positioning accuracy is improved after compensating for the errors caused by gyro-accelerometer asynchronous time. Furthermore, the gyro-accelerometer asynchronous time is caused by hardware, while the proposed method reduces the negative effect of gyro-accelerometer asynchronous time by software. Therefore, for the MIMUs with gyro-accelerometer asynchronous time, the proposed method makes it possible to apply them in PNS, which reduces the cost of pedestrian navigation and improves the reliability of MIMUs or mobile phones in PNS. Therefore, it can be concluded that the proposed calibration method works well in the foot-mounted pedestrian navigation system, and the study of gyro-accelerometer asynchronous time provides a new way to improve the positioning accuracy of pedestrian navigation in an indoor GNSS-denied environment.