1. INTRODUCTION
Pedestrian Navigation Systems (PNS) in indoor environments without GPS access can be critical in many scenarios, in particular for Search and Rescue (SAR) operations in large buildings. Certain systems rely on pre-installed infrastructures (Kuusniemi et al., Reference Kuusniemi, Bhuiyan, Strom, Soderholm, Jokitalo, Chen and Chen2012; Saad and Nakad, Reference Saad and Nakad2011; Yu et al., Reference Yu, Yang, McKelvey and Stoew2012). These can include pseudo-satellite transmitters, Radio Frequency Identification (RFID) markers and Ultra-Wideband (UWB) radars. However, the procedures of installation and maintenance of these infrastructures are too time- and energy-consuming for emergency use. Some other methods count on the geographical “fingerprint” of Radio Frequencies (RF) or geomagnetic fields (Angermann et al., Reference Angermann, Frassl, Doniec, Julian and Robertson2012; Machaj and Brida, Reference Machaj and Brida2012; Moghtadaiee and Dempster, Reference Moghtadaiee and Dempster2012), but the “fingerprints” of an area have to be acquired before positioning can be performed, which renders them inappropriate for SAR use. In addition, inertial-based systems can track the position of the pedestrian with no aid of pre-installed infrastructures and a geographical “fingerprint,” and they are generally in the form of small wearable devices, which can be worn on the helmet, belt or shoe of a person (Jimenez et al., Reference Jimenez, Seco, Prieto and Guevara2010; Patarot et al., Reference Patarot, Boukallel and Lamy-Perbal2014). Among such systems, the foot-mounted Inertial Measurement Unit (IMU) is preferred in our system, because in this configuration, the foot-mounted IMU has explicit stance phases, which can be exploited to suppress the long-term drifts of inertial sensors with the Zero-velocity UPdaTe (ZUPT)-aided Extended Kalman Filter (EKF) method (Godha and Lachapelle, Reference Godha and Lachapelle2008, Foxlin, Reference Foxlin2005; Skog et al., Reference Skog, Nilsson and Handel2010). The authors have developed a low-cost miniaturised foot-mounted pedestrian navigation system (Song et al., Reference Song, Li, Ma and Zhou2013). This also adopts a revised version of ZUPT-aided EKF, of which the threshold of zero-velocity detection is adaptively adjusted by the energy of the angular rate.
However, ZUPT-aided EKF still suffers from modelling errors, which sometimes have a more significant effect on positioning accuracy than sensor drift (Nilsson et al., Reference Nilsson, Skog and Händel2012). During the stance phases indicated by ZUPT algorithms, residual angular rate and acceleration are omitted, causing inaccurate estimation of heading and step length. During motion phases, the short-term drift of inertial sensors cannot be compensated for, and thus, it accumulates slowly. These systematic errors are intricate and are hard to eliminate in the ZUPT-aided EKF approach.
Because of the difficulty of solving the problem with sensor-oriented processing, we have to suppress the errors at a higher level as shown in Figure 1. Odometry-oriented processing only relies on “odometry,” which is constructed by positions of steps in a row. The key to odometry-oriented methods is effective exploitation of aiding knowledge on the environments to suppress error growth from a lower level ZUPT-aided EKF. Odometry-oriented methods can be categorised as the following, according to the type of knowledge that is used to aid the method:
(1) A priori knowledge of building structures. A representative method of this category is Heuristic Drift Elimination (HDE) (Borenstein and Ojeda, Reference Borenstein and Ojeda2010). HDE can correct heading drift with assumptions that there are discrete “domain” directions in most buildings, and the walking trajectories are generally aligned to these discrete directions. However, in situations with walking along complex paths, the correction procedures may result in large positioning errors, due to the wrong heading estimation, when the real heading angle significantly deviates from the “dominant” directions.
(2) Online learned location-based transition map. Pedestrian Simultaneous Localisation and Mapping (FootSLAM) creatively proposes a probabilistic transition map, which divides a regular two-dimensional (2-D) region into grids of hexagonal “cells” (Angermann and Robertson, Reference Angermann and Robertson2012). The transition map is learned online and can be exploited to aid positioning in return. FootSLAM can achieve a “bounded” positioning error with inertial sensors as sole means of navigation. However, without frequent revisiting of “cells”, FootSLAM may present particle divergence problems.
(3) Map information. A typical method exploiting map information is map-aided Cardinal Heading Aided Inertial Navigation (CHAIN) (Pinchin et al., Reference Pinchin, Hide and Moore2012). Similar to FootSLAM, this method is implemented in a particle filter framework. It uses the map of the area under investigation to “kill” impossible particles and uses the CHAIN method to achieve convergence of particles in open spaces. However, the requirements of map information are sometimes not practical in SAR applications.
Processing flow of the proposed approach.
These ingenious methods have remarkable abilities in bounding positioning errors, when they are aided by maps or characteristics of building structures. Inspired by these methods, we propose a novel approach for foot-mounted PNS that combines heading correction according to building structures with particle filter framework and introduce aiding knowledge of human behaviour to improve particle convergence. In the trajectory-calibration stage, human behaviour is learned through an inbuilt HMM with observation samples. The weight of particles can be adaptively adjusted through “discrete dominating direction” as well as the learned HMM. In that manner, a priori knowledge of building structures and human behaviour are both included in our implementation. The principle of our approach can be summarised as follows: weights of particles are determined according to the main directions of buildings and are tuned by machine learning algorithms based on human walking behaviour. The overall flow of our approach is shown in Figure 1. As the ZUPT-aided EKF in the trajectory-generation stage is already commonly used in inertial-based PNS, it is not elaborated on in this paper.
The real trajectory of a pedestrian is three dimensional (3-D), but the vertical position resolved by ZUPT-aided EKF may be inaccurate. The authors have already proposed an approach for accurate height estimation (Gu et al., Reference Gu, Ma, Li and Song2013), which makes use of a priori knowledge of building structures to reduce the drifts of inertial sensors and barometers. As a result, this paper only investigates positioning in planar dimensions.
2. THE PARTICLE FILTER FRAMEWORK
2.1. Particle filter algorithm
In contrast to Gaussian Filters, the particle filter approach approximates the real posterior density by finite samples, rather than a fixed function form. This makes the particle filter more suitable for complex posterior representations without making too many assumptions on the function's parameters (Doucet and Johansen, Reference Doucet and Johansen2009). The particle filter framework is adopted in our implementation for two reasons. (1) The Probability Density Function PDF) of positions are non-linear and non-stationary. (2) It is flexible in exploiting aiding knowledge. These reasons are elaborated separately in detail as follows.
As an odometry-oriented method, the particle filter only deals with position vectors from lower level results. Coupled gyroscope errors and accelerometer errors are propagated in a rotating reference frame, causing complicated time-varying probabilistic distribution of positions (Titterton and Weston, Reference Titterton and Weston2004. P.342). Under such circumstances, particle filters, as a nonparametric implementation, are suitable for representing the time-varying distribution. Another reason for choosing the particle filter approach is that the human motion model is highly non-linear. The non-linear versions of the Kalman Filter, such as EKF and Unscented Kalman Filter (UKF), handle non-linear models by linearizing state transitions through a Taylor expansion or calculating the linearized approximation at selected points. However, the linearization may be a poor approximation of the state transition function, which induces additional errors. In that manner, a particle filter performs better by directly calculating state samples or particles, which matches our situation better.
In addition to the inherent advantages of a particle filter, it fits our application for another reason: the flexibility in incorporating aiding knowledge. In the particle filter framework, particles can be “killed,” de-weighted or over-weighted according to aiding knowledge. For example, when a particle appears in a wall, it is “killed” because the map-conditioned probability is zero at that spot.
Particle filter methods have been widely adopted recently in the field of indoor navigation. However, the way in which aiding knowledge is used at the weight update stage differs in various methods. Commonly used aiding knowledge includes map information or a priori knowledge of building structures, but human behaviour is often omitted. In this paper, both a priori knowledge of building structures and human walking behaviour are integrated into the proposed approach to update weight of particles, and there are significant performance effects in the trajectory-calibration stage.
2.2. Particle filter stages
2.2.1. Definitions
Odometry measurements are used to describe the movement of a subject over time. Unlike classical robotic motion models (Thrun et al., Reference Thrun, Burgard and Fox2005. P.107), odometry measurements in human motion models are embedded in a fixed navigation reference frame and not in the internal reference frame of the subject. The reasons for exploiting fixed navigation reference frame embedded measurements as inputs to the particle filter are simple. (1) The results of lower-level ZUPT-aided EKF resolution are already navigation reference frame embedded. (2) A pedestrian's gaits are sometimes irregular, including movements such as side-stepping, where the pedestrian's facing (the direction of an object's front) deviates greatly from the heading (the direction in which an object moves) during a step. This makes it harder for the pedestrian's internal reference frame to reflect motions than a fixed one. (3) Aiding knowledge can be exploited more conveniently in a fixed reference frame. During the stance phase in each step, the step position (x t,yt) is acquired from the trajectory-generation stage. The odometry measurement, which contains a translation and a rotation as shown in Figure 2, can be calculated from the step position as follows.
$$\eqalign{\Delta x_t & = x_t - x_{t - 1} \cr \Delta y_t & = y_t - y_{t - 1} \cr \theta _t & = \arctan \left( {\displaystyle{{\Delta y_t} \over {\Delta x_t}}} \right) \pm \pi \cr \delta _{rotation,t} & = \Delta \theta _t = \theta _t - \theta _{t - 1} \cr \delta _{translation,t} & = l_t = \sqrt {\Delta x_t ^2 + \Delta y_t ^2}} $$Human motion model.
Though variable set {x t,yt,Δx t,Δy t,θt,Δθ t,lt} has many compact subsets, we choose measurement Z t = {δ rotation,t,δtranslation,t} = {Δθ t,lt} as an input to the particle filter for the purpose of easy implementation.
Similar to measurements, the state space of a pedestrian may only be comprised of two components x r,t and y r,t. However, we can add redundancies θ r,t and l r0 to make the recursion easier, and the suffix r indicates the value of variable is real.
To start the state recursion, the initial state should be decided. It is closely connected to the choice of navigation reference frame. To make it simple, we choose the IMU initial position to be the origin of the navigation reference frame and the IMU initial facing to be the x-axis of it. When the IMU initial position is previously acquired through GPS, and when the angle between the IMU initial heading and north is acquired through a compass, the navigation frame can be easily transformed to the Earth frame by translation and rotation. As a result, the initial state can be denoted as follows:
$$ S_0 = \left\{ {\matrix{ {x_{r0}} \cr {y_{r0}} \cr {\theta _{r0}} \cr {l_{r0}} \cr}} \right\} = \left\{ {\matrix{ 0 \cr 0 \cr 0 \cr 0 \cr}} \right\}$$2.2.2. Particle propagation
The new state particle set at time t should be proportional to the posterior S t[1:N] ~ p(S t|S t−1,Z t), where N denotes the total number of particles. We assume that the noise underlying p(S t|St− 1,Z t) is made up of two independent white Gaussian processes that can be represented in Equation (3).
$$\eqalign{\,p(l_{r,t} \vert l_t )& =\displaystyle{1 \over {2\pi \sigma _l}} \exp \left( { - \displaystyle{{(l_t - l_{r,t} )^2} \over {2\pi \sigma _l ^2}}} \right) \cr p(\Delta \theta _{r,t} \vert \Delta \theta _t ) &=\displaystyle{1 \over {2\pi \sigma _{\Delta \theta}}} \exp \left( { - \displaystyle{{(\Delta \theta _t - \Delta \theta_{r,t} )^2} \over {2\pi \sigma_{\Delta \theta}^{2}}}} \right)} $$where σ l2 denotes the variance of translation noise, and σ Δθ2 denotes the variance of rotation noise. To make the noise easier to analyse, we assume that σ l2 and σ Δθ2 are proportional to the absolute value of l t and ∆θ t in Equation (4).
$$\eqalign{\sigma _l ^2 &= c\left\vert {l_t} \right\vert \cr \sigma _{\Delta \theta} ^2 &= d\left\vert {\Delta \theta _t} \right\vert } $$where c and d are scale coefficients. The model in Equation (4) may not be accurate, but from our observations, it can reveal the phenomenon that greater translation and rotation generally lead to larger noise. In this case, the posterior function p(S t|St− 1,Z t) is hard to represent as an analytic formulation. As a result, the particle set at time t is not directly sampled from p(S t|St− 1,Z t) but is derived in Table 1.
Algorithm for particle derivation with time.
2.2.3. Weight Update and Resampling
In the conventional particle propagation process, the particles are all equally weighted. Without any weight update or resampling process applied, the state particle set S t[1:N] would diverge quickly, resulting in poor algorithm robustness. We implement a novel weight update strategy in our approach, where the weight of particles is adaptively adjusted according to a priori knowledge of building structures and human behaviours. We only illustrate how building structure is exploited in this section, whereas human behaviour learning will be depicted individually in the next section.
In most buildings, the corridors are straight and are either parallel or orthogonal to each other. Even in some rooms, the layouts are consistent with the four “dominant” directions of corridors. This common characteristic of buildings can be exploited to improve particle convergence. As shown in Figure 3 (a), eight discrete directions are predefined according to the layout of the building. d1, d3, d5 and d7 are the directions of corridors in the building, whereas d2, d4, d6 and d8 are complementary directions for the situation of complex trajectories, for example, walking causally in a room. The complementary directions make the direction adjustment less arbitrary. Thus, particles near the eight directions are assigned weights according to Gaussian functions as follows. We call this weight assignment rule the Discrete Dominant Direction (DDD) rule in this paper. ψ 1:8 are the angles of the eight directions in the navigation reference frame.
(a) Heuristic directions in buildings. (b) Weight update strategy.
Algorithm for particle weight updating, assuming rectangular structure of buildings.
In the algorithm, a 1 to a 8 are coefficients of the Gaussian functions and can be considered as the “weights” of the eight directions as shown in Figure 3(b). They are determined by human behaviour learning and will be explained in the next section.
The resampling algorithm we use in our implementation is residual resampling (Hol et al., Reference Hol, Schon and Gustafsson2006). It is used to eliminate particles with small importance weights and duplicate particles with significant weights.
3. HUMAN WALKING BEHAVIOR LEARNING BASED ON HMM
Here, human walking behaviour can be uncovered by the pedestrian's odometry measurements. There are many limitations on the trajectory according to the pedestrian's walking behaviours. For example, it often takes 3~5 steps to complete a 90° or 180° turn. Straight lines are connected with bent turning curves. Translation during a step is generally shorter when the pedestrian is about to take a turn or during a turn. However, information involving human walking behaviour is not effectively exploited in any of the existing methods, which are all under the Markov assumption that the next step is only relative to the current one but irrelevant to previous ones, which interrupts connections between several neighbouring steps. To resolve the problem, a modified HMM algorithm is proposed here to extract human walking behaviours. To make our method more applicable, we can consider that human walking behaviour is common in nature among people, for example, a pedestrian can walk straight ahead or turn left/right ahead, but he or she would rarely walk straight backward or turn left/right backward. Similarly, the structure of the building is not limited to a single given building, because this regular structure, such as right angles and long straight lines, is generally seen in almost all buildings.
3.1. HMM
The key feature in a pedestrian's trajectory is measured rotation Δθ t in each step. Δθ t is sorted into eight discrete states values according to the region it falls into as shown in Figure 4 (note that φ 1 to φ 8 represent the rotation angle set {0, π /4, π/2, 3π/4, π, − 3π/4, −π/2, −π/4}). The difference between Figure 3(a) and Figure 4 is that in Figure 3(a), the angles ψ 1:8 represent angles of eight fixed directions in the navigation reference frame, but, in Figure 4, φ 1:8 represent angles of eight directions in rotation Δθ t, φ 1:8 change over time in the navigation reference frame.
Participation of Δθ t to eight regions.
We define a function r s=Region(θ), where the input θ is an angle, and output r s is the region θ falls into. If we only choose a single Region(Δθ t) to represent a state in HMM, the transition can only manifest a relationship between two adjacent steps. Here, we put together three of these rotations to form a new state.
The new states Y t are taken from a sliding window and overlap each other as in Figure 5. Although Y t is able to represent a sub-trajectory composed of three steps, it causes some issues because the number of its possible values can reach 8 × 8 × 8=256. The great number of results in the state transition matrix and observation matrix are large dimensioned, which makes it harder for further training to achieve convergence. To solve this problem, simple cluster analysis was performed to lower the number of state values. According to our observation comprising 50 trajectories, some of the state values never emerge, and some of the values emerge at a very low frequency. Those never-shown values are eliminated from possible state values, whereas those seldom shown (lower than 1%) values are sorted out into the “abnormal” category with the value 0. With the clustering, the amount of state values decrease to only 14. They are {r 1r 1r 1,r 1r 1r 2,r 2r 1r 1,r 8r 1r 1,r 1r 1r 8,r 1r 2r 2,r 2r 2r 1,r 2r 8r 1,r 1r 8r 8,r 8r 8r 1,r 1r 2r 8,r 1r 8r 1,r 1r 2r 3,0}. To simplify the value, they are replaced by {q 1,q 2,…,q 14}, respectively.
State definition in HMM model
The states previously discussed are observations, whereas the hidden states X t are defined exploiting real rotation Δθ r,t, which has a one-to-one mapping relationship to observations.
3.2. Training parameters of HMM
After the foregoing sequences of observation are acquired, we can estimate the parameter λ of HMM from a learning process. After learning, the new parameter λ new would explain the observation sequences at least as well as or better than the initial estimation λ initial, such that p(Y|λ new)≥p(Y|λ initial). To perform the unsupervised learning, we use the Baum-Welch algorithm (Welch, Reference Welch2003; Poritz, Reference Poritz1988), which make use of an Expectation Maximisation (EM) algorithm to find the maximum likelihood estimation of. λ.
The value set of states and observation is {q 1,q 2,…,q 14} and the history sequence of observation is {o 1,o 2,…,o t,ot +1,….} The parameter to be trained is λ=(A,B,π), where A and B denote the transition probability matrix and the emission probability matrix, respectively, and π denotes the initial probability array as shown in Equation (5).
$$\eqalign{& A = [a_{ij} ] = p(X_t = q_i \vert X_{t - 1} = q_j ) \cr & B = [b_{ij} ] = p(Y_t = q_j \vert X_t = q_i ) \cr & \pi _i = p(X_1 = q_i )} $$The algorithm has three procedures:
(1) Forward procedure: The forward variable is α i(t)=p(Y 1=o 1,…,Y t=o t,Xt=q i|λ), which denotes the given HMM parameters, the probability of being at the state q i and having observed sequence {o 1,o 2,…,o t} at time t. b i(o t) denotes p(Y t=o t|Xt=q i). α i(t). is recursively derived in Equation (6).
(6)
$$\eqalign{\alpha _i (1) &= \pi _i b_i (o_1 ) \cr \alpha _j (t + 1) &= b_j (o_{t + 1} )\sum\limits_{i = 1}^N {\alpha _i (t)a_{ij}}} $$(2) Backward procedure: The backward variable is β i(t)=p(Y t+1=o t+1,…,Y T=o T|Xt=q i,λ), which denotes the given HMM parameters and the fact of being in state q i at time t, the probability of going to observe a sequence {o t+1,o t+2,…,o T}. βi(t) is recursively derived in Equation (7).
(7)$$\eqalign{\beta _i (T) &=1 \cr \beta _i (t) &=\sum\limits_{\,j = 1}^N {\beta _j (t + 1)a_{ij} b_j (o_{t + 1} )}} $$(3) Update procedure: Temporary variables are introduced in Equation (8).
(8)$$\eqalign{& \gamma _i (t) = p(X_t = q_i \vert Y,\lambda ) = \displaystyle{{\alpha _i (t)\beta _i (t)} \over {\sum\limits_{\,j = 1}^N {\alpha _j (t)\beta _j (t)}}} \cr & \xi _{ij} (t) = p(X_t = q_i, X_{t + 1} = q_j \vert Y,\lambda ) = \displaystyle{{\alpha _i (t)a_{ij} \beta _j (t + 1)b_j (o_{t + 1} )} \over {\sum\limits_{k = 1}^N {\sum\limits_{l = 1}^N {\alpha _k (t)a_{kl} \beta _l (t + 1)b_l (o_{t + 1} )}}}}} $$
$\sum\nolimits_{t = 1}^{T - 1} {\gamma _i (t)} $ can be interpreted as the expected number of transitions from q i, and 
$\sum\nolimits_{t = 1}^{T - 1} {\xi _{ij} (t)} $ can be interpreted as the expected number of transitions from q i to q j. Based on the expectations, the parameters are updated in the following steps.
(1) πi′=γ i (1) (expected frequency for state q i to appear at time t=1)
(2)
$a_{ij} \hskip-3pt\prime = {{\sum\nolimits_{t = 1}^{T - 1} {\xi _{ij} (t)}} / {\sum\nolimits_{t = 1}^{T - 1} {\gamma _i (t)}}} $ (expected number of transitions from q i to q j divided by expected number of transitions from q i)(3)
$b_{ij}\hskip-3pt\prime = {{\sum\nolimits_{t = 1,{\rm } s.t.{\rm } o_t = q_j} ^{T - 1} {\gamma _i (t)}} / {\sum\nolimits_{t = 1}^{T - 1} {\gamma _i (t)}}} $ (expected number of times in state q i and observing q j divided by expected number of times in state q i)
The procedures are iterated until reaching convergence. Because the hidden states have a one-to-one mapping relationship with observation, the initial guess of A is in Equation (9).
$$a_{ij} ^0 = \displaystyle{{number\, of\, transitions\, from\, q_j \, to\, q_i \, in\, observation\, sequence} \over {number\, of\, times\, in\, state\, q_i \, in\, observation\, sequence}}$$The initial guess of B is based on our experience that rotation ∆θ t has a relatively small error in each step, so 
$b_{ij}^0 $ is set to 0.9, and the left over 0.1 is equally allocated to other “adjacent” b ik. (b ij and b ik are “adjacent” means q j and q k are “adjacent” as shown in Figure 6).
Definition of “adjacent” regions.
3.3. Weight update strategy based on parameters of HMM
Once the parameters are decided, we use A, B and the previous observation Y t−1 to predict the next state X t. If Y t−1=q 1, then the prediction can be denoted as:
$$\eqalign{ p(X_t = q_k \vert Y_{t - 1} = q_l ) &= \sum\limits_{\,j = 1}^N {\,p(Y_{t - 1} = q_l \vert X_{t - 1} = q_j )p(X_t = q_k \vert X_{t - 1} = q_j )} \cr & = \sum\limits_{\,j = 1}^N {b_{\,jl} a_{kj}}} $$p(X t=q k|Yt− 1=q l) is the probability of the state X t equals one of the 14 states values {q 1,q 2,…,q 14}, where the states are composed of three steps. However, the weight update process is aimed at update w t[1:N] in a single step, so p(X t=q k|Yt− 1=q l) have to be converted to p(Region(Δθ t)=r i|Yt− 1=q l), which denotes the probability that Δθ t falls into each region r i. The conversion is conducted in Equation (11). (The “abnormal” state is equally allocated to r i.)
$$\eqalign{& p({\rm Region(}\Delta \theta _t {\rm )} \!=\! r_1 {\rm \vert }Y_{t - 1}\!=\! q_l {\rm )}\!=\! \sum\limits_{k = 1,3,4,7,8,10,12} {\,p(X_t \!=\! q_k \vert Y_{t - 1} \!=\!q_l )}\cr&\hskip138pt+ p(X_t \!=\! q_{14} \vert Y_{t - 1} \!=\! q_l )/8 \cr & p({\rm Region(}\Delta \theta _t {\rm )} \!=\! r_2 {\rm \vert }Y_{t - 1} \!=\! q_l {\rm )} \!=\! \sum\limits_{k = 2,6} {\,p(X_t \!=\! q_k \vert Y_{t - 1} \!=\! q_l )} + p(X_t \!=\! q_{14} \vert Y_{t - 1} \!=\! q_l )/8\cr & p({\rm Region(}\Delta \theta _t {\rm )} \!=\! r_3 {\rm \vert }Y_{t - 1} \!=\! q_l {\rm )} \!=\! p(X_t \!=\! q_{13} \vert Y_{t - 1} \!=\! q_l ) + p(X_t \!=\! q_{14} \vert Y_{t - 1} \!=\! q_l )/8 \cr & p({\rm Region(}\Delta \theta _t {\rm )} \!=\! r_8 {\rm \vert }Y_{t - 1} \!=\! q_l {\rm )} \!=\!\sum\limits_{k = 5,9,11} {\,p(X_t \!=\! q_k \vert Y_{t - 1}\!=\! q_l )} + p(X_t \!=\! q_{14} \vert Y_{t - 1} \!=\! q_l )/8 \cr & p({\rm Region(}\Delta \theta _t {\rm )} \!=\! r_{4,5,6,7} {\rm \vert }Y_{t - 1} \!=\! q_l {\rm )} \!=\! p(X_t \!=\! q_{14} \vert Y_{t - 1} \!=\! q_l )/8} $$The equations denote that the relations between the state values {q 1,q 2,…,q 14} and the regions {r 1,r 2,…,r 8}, for example, q 2={r 1,r 1,r 2} and q 6={r 1,r 2,r 2}, they are the only two of the 13 state values (except q 14, the “abnormal” state is equally allocated to r i) which end up with region r 2, so the equation
$$\eqalign{ p({\rm Region(}\Delta \theta _t {\rm )} = r_2 {\rm \vert }Y_{t - 1} = q_l {\rm )} &= \sum\limits_{k = 2,6} {\,p(X_t =q_k \vert Y_{t - 1} = q_l )} \cr&\quad+ p(X_t = q_{14} \vert Y_{t - 1} = q_l )/8\ {\rm is\, true}$$Assuming that p s=p(Region(Δθ t)=r s|Yt− 1=q l), then we can obtain the predicted region probability set {p 1,p 2,…p 8}.
The region probability set is used to determine the coefficients a 1 to a 8 of the weight update algorithm described in the second section. As shown in Figure 7, each region r i has a corresponding direction φ i. The region probability set {p 1,p 2,…p 8} can be counted as the probabilities of the discrete directions{φ 1, φ 2,…, φ 8}. The probabilities of these directions are projected to the eight fixed heuristic directions to determine the coefficients a 1 to a 8.
An example of projections from predicted directions to heuristic directions.
Given the heuristic directions and the predicted directions of Δθ t, the coefficients are the sum of projections of adjacent φ. As shown in the example in Figure 7, we can see that a 1=p 1 · cos(π/4−α) + p 8 · cos(α). After the coefficients are decided, the weights of particles are assigned according to Figure 3(b).
4. TESTS AND EXPERIMENTS
Based on the miniaturised foot-mounted pedestrian navigation system (Song et al., Reference Song, Li, Ma and Zhou2013), the training samples were acquired by two people walking in a typical office building. The 15 walks cover different scenarios, for example, corridors, rooms with restricted space, rooms with large open space and stairs. It can be considered that the training is basically independent from the subject/environment and that the trained parameters could be used in other environments. However, to acquire more accurate model parameters, we still need more samples in the training process. In our future work, we will collect samples from more people in different buildings, and investigate how training samples affect the overall performance, but now we just exploit the samples we considered adequate to perform the training offline. After the training, trained parameters of HMM are exported to the implementation algorithm. In all the following experiments, the number of particles is chosen to be 1000. This number is chosen because the offline processing time is acceptable and the location accuracy would not improve much by further increasing the particle number. It is a compromise between time consumed and location accuracy.
To verify the accuracy and robustness of our approach, we conducted a large-scale outdoor experiment in an extreme situation, including many sharp turns, curves and straight line trajectories, and it is worth noting that there are few location revisits, which means no foregoing information can be adopted for trajectory calibration. The start and end point of this route are the same, so the final error of positioning can be easily calculated. The duration of the route is approximately 12 minutes, and the total distance is approximately 908 metres in 699 steps (sensor on one foot only). Without any aiding maps or location revisits, our approach achieves an error of no more than 2.5 metres, whereas the pure ZUPT-aided EKF and HDE method shows much larger errors, as shown in Figure 8. The pure ZUPT-aided EKF method has significant heading drifts, which severely deteriorate location accuracy. The HDE method can suppress heading drifts in certain directions, but when the trajectories are curves, it may lead to divergence. The FootSLAM method is not compared here, because that, by its principle, requires frequent location revisits, as otherwise it may degenerate into a normal particle filter method. Map-aided CHAIN is not compared here either because the map information is not known in our application.
Trajectory comparison of three methods in an outdoor walking experiment.
Another experiment is conducted in a typical office building, which contains real indoor scenarios, including rooms, corridors and stairs. The people walk upstairs to the sixth floor and then go downstairs to the start point, during which several corridors have been crossed as shown in Figure 9. The 3-D and 2-D trajectories reveal good consistency in our approach. In the projected 2-D plane, the trajectories in different floors coincide well with each other. In particular, in each of the three stair regions (red circles), the error maintains a range of less than 1 metre. This result has validated the effectiveness of our method in bounding errors during the walk. Our approach only copes with 2-D plane positions because the accurate estimation of height has already been resolved by the authors, based on the intelligent recognition of stairs and corridors (Gu et al., Reference Gu, Ma, Li and Song2013).
Indoor multi-floor consistency test with rooms, corridors and stairs.
Like any inertial methods based on dead reckoning principles, the current positioning accuracy of foot-mounted PNS is highly dependent on the accuracies of previous locations and headings. If fairly good accuracy is not maintained after a small-scale complex walk, long-term robustness also suffers. Here, we conduct an experiment to evaluate the performance under a complex trajectory situation in Figure 10. During our test, the trajectory is highly irregular with frequent random curves from room to room. As observed in the results, the last segment of the trajectory is parallel to the first segment, which denotes the right estimation of direction. In addition, the total error after the test is less than 1 metre, which is acceptable in most SAR uses.
Complex trajectory with frequent direction changes.
The authors have recently conducted another experiment to test the performance of our method when the pedestrian makes unexpected movements. In the experiment, the pedestrian walks through a line segment back and forth. At the endpoints of the line segment, the pedestrian moves backwards suddenly, which can be regarded as the “abnormal” state q 14. As shown in the trajectories from Figure 11, the trajectory of our method almost coincides with that of pure ZUPT-aided EKF. Although the error of our method, 0.92 m, is larger than pure ZUPT-aided EKF 0.5 m, it is still acceptable considering unexpected motions are rare in a pedestrian's daily activities.
Comparisons of the methods when a pedestrian makes unexpected movements.
It is hard to design a method to find a fixed relationship between the total walking distance and error growth, especially when the trajectories are too complex. However, we can dissect the error growth in simpler situations, when the trajectory is in simple forms, for example, a square in Figure 12. We conducted an experiment along a square with 4 m side, and the results indicate that our approach has good performance in bounding errors of both location and heading.
Error growth comparison when walking around a square.
5. CONCLUSIONS
In this paper, a novel approach for foot-mounted inertial PNS is proposed with few presuppositions and requirements. For the method, a priori knowledge on human walking behaviour and building structures are innovatively integrated in our approach to update the weight of particles adaptively. The key to our approach lies in the weight update strategy. Particles are weighted according to the DDD rule and the determination of particle weight is through an HMM whose parameters are learned from human walking behaviour.
The results of experiments on different indoor scenarios have demonstrated the superiority of our approach. Its advantages can be summarised as follows: (1) good accuracy and consistency when walking simple trajectories; (2) location revisits are not a prerequisite; (3) adaptability to some special situations with complex paths and unexpected movements.
The integration of the DDD rules and human behaviour achieved better results than pure ZUPT-aided EKF and HDE. In addition, compared with FootSLAM and map-aided CHAIN, our method can be applied to harsh applications, such as SAR operations, very well, where frequent location revisits and map information are often not available.
