On Monte Carlo Smoothing in Multi Sensor Indoor Localisation

Toni Fetzer, Frank Ebner, and Frank Deinzer

Faculty of Computer Science and Business Information Systems

University of Applied Sciences Würzburg-Schweinfurt

Würzburg, Germany

{toni.fetzer, frank.ebner, frank.deinzer}@fhws.de

Lukas Köping, Marcin Grzegorzek

Pattern Recognition Group University of Siegen Siegen, Germany {lukas.koeping, marcin.grzegorzek}@uni-siegen.de

Abstract—Indoor localisation continuous to be a topic of growing importance. Despite the advances made, several profound problems are still present. For example, estimating an accurate position from a multimodal distribution or recovering from the influence of faulty measurements. Within this work, we try to solve such problems with help of Monte Carlo smoothing methods, namely forward-backward smoother and backward simulation. In contrast to normal filtering procedures like particle filtering, smoothing methods are able to incorporate future measurements instead of just using current and past data. This enables many possibilities for further improving the position estimation. Both smoothing techniques are deployed as fixed-lag and fixed-interval smoother and a novel approach for incorporating them easily within our localisation system is presented. All this is evaluated on four floors within our faculty building. The results show that smoothing methods offer a great tool for improving the localisation results. Especially fixed-lag smoothing provides a great runtime support by reducing timely errors and improving the overall estimation with affordable costs.

I. Introduction

Determining a position indoors is a challenging task. Besides the complex architecture of many buildings, high accuracy needs to be achieved, especially for buildings with many small separated areas like shopping malls or office blocks. In recent years, many different systems were presented to meet those requirements. Especially Wi-Fi positioning and pedestrian dead reckoning (PDR) are very popular solutions. Approaches based on PDR try to estimate the current position given the previous position and thus require an initial state. However, this allows for cumulative errors and leads to an erroneous position estimation within a very short period. By incorporating the absolute position information of Wi-Fi this drift can be corrected. Additional improvements can be achieved by using environmental information about walls and obstacles provided by a floor map.

In most cases, probabilistic methods are used to incorporate those highly different sensor types. Here, a probability distribution describes the pedestrian’s possible whereabouts and therefore the uncertainty of the system. Drawing from a probability distribution and finding an analytical solution for densities is in most cases a difficult task, especially in case of time sequential, non-linear and non-Gaussian models. Due to the high complexity of the human movement, we consider indoor localisation as such.

Bayesian filters solve such problems by updating the state estimation recursively with every new incoming measurement. A broad class to obtain numerical results for this approach are Monte Carlo (MC) methods, where a set of random samples is used to approximate the underlying probability distribution. By applying the time sequential hidden Markov process of Bayes filtering, one of the most important MC techniques results: particle filtering. Here, a set of weighted random samples is used to solve the state estimation process.

Based on this general methodology, many different approaches for estimating a position in indoor environments have been developed. All these approaches differ mainly in how the dynamics are modelled in the transition step and how a specific sensor measurement can be used for evaluation. For example, recent approaches are using a graph-based structure to consider environmental restrictions (walking through walls) and the characteristics of human movement (walking speed) within the transition model [1]–[3]. The evaluation model is mostly separated into any number of sensor models, each representing the probability for a noisy measurement in regard to the current position. For example, a barometer can be used to determine the probability of being on a certain floor [4].

Despite the many advances made in the last years, nearly all systems suffer from more or less the same problems. Like mentioned before, PDR suffers from an accumulating bias, the signal of Wi-Fi gets attenuated by walls and the barometric pressure is highly affected by weather patterns and humidity [4]. That is the reason for the use of statistical methods in the first place. Nevertheless, there are even more profound problems regarding the whole position estimation procedure.

Current transition models, which aim to approximate the movement, are still very restrictive and unable to handle unforeseen events. Faulty sensor measurements, like a falsely detected turn, can cause the estimation to lose track. For example by recognising a turn too soon and walking into a room instead of another big hallway. Due to this, the filter needs some time to recover, which again takes a while because of the restrictive model (e.g. no walking through walls and only realistic walking speed). This temporal delay worsens the estimate immensely. A solution to recover from such

Time t  t − 1  t − 2

Fig. 1. An example of the occurrence of a multimodal distribution. At time t − 1 the floor is separated by a wall and the distribution (coloured circle), splits apart. The most likely position (green line) is estimated somewhere in-between. After a right turn, the distribution slowly starts to recover its unimodality.

filter divergences faster, involves methods to re-initialize the filtering procedure [2]. However, even this can not completely prevent delays. Another reason for possible time delays are slow sensor updates. For example, most mobile devices restrict the Wi-Fi module to update only every few seconds, to save on battery. Further critical problems arise from multimodal distributions. Those are caused by multiple possible position estimates. Fig. 1 illustrates an example where a floor gets separated by a wall. Due to inaccurate measurements and a PDR approach for evaluating the movement, the distribution splits apart. Therefore, the weighted average position is somewhere in-between. Only after the pedestrian turns right, the distribution is again unimodal, since moving through walls is impossible. As one can imagine, this can lead to serious problems in big indoor environments. Such a situation can be improved by incorporating future measurements (e.g. the right turn) or predictive information (e.g. the most likely path) to the filtering procedure [5]. However, standard filtering methods are not able to use any future information and the possibilities to make a distant forecast are also limited [6]–[8].

One very promising way to deal with these problems is smoothing. Smoothing methods are able to make use of future measurements for computing their estimation. By running backwards in time, they are also able to remove multimodalities and improve the overall localization result. Since the problem of navigation, especially the representation of complex movement patterns, results in a non-linear and non-Gaussian state space, this work focuses mainly on smoothing techniques based on the broad class of MC methods. Namely, forward-backward smoothing [9] and backward simulation [10].

Within this work, we investigate the benefits and drawbacks of those techniques using our indoor localisation system presented in [5]. We provide both, fixed-lag and fixed-interval smoothing as well as a novel approach for incorporating them easily within the localisation procedure. Additionally, we enrich the state transition model with an activity recognition to distinguish between walking, standing and walking stairs. The main goal is to solve above mentioned problems and to investigate new possibilities for even more advanced systems.

All of our contributions are supported by an extensive experimental evaluation.

II. Related Work

Sequential MC filter, like aforementioned particle filter, use all observations o_{1:t} until the current time t for computing an estimation of the state q_{t}. In a Bayesian setting, this can be formalized as the computation of the posterior distribution p(q_{t} | o_{1:t}) using a sample of N independent random variables, X_{t}^{i} ∼ (q_{t} | o_{1:t}) for i = 1, ..., N for approximation. Due to importance sampling, a weight W_{t}^{i} is assigned to each sample X_{t}^{i}. In context of particle filtering {W_{1:t}^{i}, X_{1:t}^{i}}_{i=1}^{N} is a weighted set of samples, also called particles. Therefore a particle is a representation of one possible system state q. By considering a situation given all observations o_{1:T} until a time step T, where t ≪ T, standard filtering methods are not able to make use of this additional data for computing p(q_{t} | o_{1:T}). This problem can be solved with a smoothing algorithm.

Within this work we utilise two types of smoothing: fixedlag and fixed-interval smoothing. In fixed-lag smoothing, one tries to estimate the current state, given measurements up to a time t+τ, where τ is a predefined lag. This makes the fixed-lag smoother able to run online. On the other hand, fixed-interval smoothing requires all observations until time T and therefore only runs offline, after the filtering procedure is finished [8].

The origin of MC smoothing can be traced back to Genshiro Kitagawa. In his work [11] he presented the simplest form of smoothing as an extension to the particle filter. This algorithm is often called the filter-smoother since it runs online and a smoothing is provided while filtering. This approach uses the particle filter steps to update weighted paths {(W_{t}^{i}, X_{1:t}^{i})}_{i=1}^{N}, producing an accurate approximation of the filtering posterior p(q_{t} | o_{1:t}) with a computational complexity of only O(N). However, it gives a poor representation of previous states due a monotonic decrease of distinct particles caused by resampling of each weighted path [7]. Based on this, more advanced methods like the forward-backward smoother [9] and backward simulation [10] were developed. Both methods are running backwards in time to reweight a set of particles recursively by using future observations. Algorithmic details will be shown in section IV.

In recent years, smoothing got attention mainly in other areas than indoor localisation. The early work of [12] demonstrates the possibilities of smoothing for visual tracking. They used a combination of the CONDENSATION particle filter with a forward-backward smoother. Based on this pioneering approach, many different solutions for visual and multi-target tracking have been developed [13]. For example, in [14] a particle smoother is used to reduce multimodalities in a blood flow simulation for human vessels. The authors of [15] use a smoother to overcoming the problem of particle impoverishment while predicting the Remaining Useful Life (RUL) of equipment (e.g. a Lithium-ion battery).

Nevertheless, there are some promising approaches for indoor localisation systems as well. For example [2] deployed a fixed-interval forward-backward smoother to improve the position estimation for non-real-time applications. They combined Wi-Fi, step and turn detection, a simple line-of-sight model for floor plan restrictions and the barometric change within a particle filter. The state transition samples a new state based on the heading change, altitude change and a fixed step length. The experiments of [2] clearly emphasize the benefits of smoothing techniques. The estimation error could be decreased significantly. However, a fixed-lag smoother was discussed only in theory.

In the work of [16] both fixed-interval and fixed-lag smoothing were presented. They implemented Wi-Fi, binary infra-red motion sensors, binary foot-switches and a potential field for floor plan restrictions. Those sensors were incorporated using a sigma-point Kalman filter in combination with a forwardbackward smoother. It was also proven by [16], that the fixedlag smoother is slightly less accurate than the fixed-interval smoother, as one would expect from the theoretical foundation. Unfortunately, even a sigma-point Kalman filters is after all just a linearisation and therefore not as flexible and suited for the complex problem of indoor localisation as a non-linear estimator like a particle filter. Additionally, the Wi-Fi RSSI model requires known calibration points and is deployed using a remarkable number of access points for very small spaces. In our opinion this is not practical and does not suite realworld conditions. Since humans with a specific destination in mind do not tend to change their directions randomly, we would further recommend adding a PDR-based transition to draw samples in a more directed manner instead of scattering them randomly in every direction.

The here presented approach will use two different smoothing algorithm, both implemented as fixed-interval and fixed-lag versions. Further, our localisation system presented in [5] enables us to walk stairs and thus going into the third dimension. Therefore, a regularly tessellated graph is utilized to avoid walls, detecting doors and recognizing stairs. Within this work, this is additionally supported by a simple classification that detects the activities unknown, standing, walking and walking stairs. Finally, we incorporate prior navigation knowledge by using syntactically calculated realistic human walking paths [5]. This method makes use of the given destination and thereby provides a more targeted movement.

III. Filtering

As mentioned before, most smoothing methods require a preceding filtering. In our previous work [5], we consider indoor localisation as a time-sequential, non-linear and non-Gaussian state estimation problem. Therefore, a Bayes filter that satisfies the Markov property is used to calculate the posterior, which is given by

p(q_{t} | o_{1:t}) ∝

\underbrace{p(o_{t} | q_{t})}_{evaluation}∫ \underbrace{p(q_{t} | q_{t−1}, o_{t−1})}_{transition} \underbrace{p(q_{t−1} | o_{1:t−1})dq_{t−1}}_{recursion} .

(1)

Here, the previous observation o_{t−1} is included into the state transition [17]. For approximating eq. (1) by means of MC methods, the transition is used as proposal distribution, also known as CONDENSATION algorithm [12]. The handle the phenomenon of weight degeneracy we additionally apply a resampling step.

In context of indoor localisation, the hidden state q is defined as follows:

q = (x, y, z, θ, ρ̂_{rel}), x, y, z, θ, ρ̂_{rel} ∈ ℝ ,  (2)

where x, y, z represent the position in 3D space, θ the user’s heading and ρ̂_{rel} the relative atmospheric pressure prediction in hectopascal (hPa). Further, the observation is given by

o = (s_{wifi}, s_{ib}, Δθ, n_{steps}, ρ_{rel}, Ω) ,  (3)

covering all relevant sensor measurements. Here, s_{wifi} and s_{ib} contain the measurements of all nearby access-points (AP) and iBeacons, respectively. Δθ and n_{steps} describe the relative angular change and the number of steps detected for the pedestrian. ρ_{rel} is the relative barometric pressure with respect to a fixed reference. Finally, Ω contains the activity, currently estimated for the pedestrian, which is one of: unknown, standing, walking, walking stairs up or walking stairs down. The probability density of the state evaluation is given by

p(o_{t} | q_{t}) = p(o_{t} | q_{t})_{baro} p(o_{t} | q_{t})_{ib} p(o_{t} | q_{t})_{wifi} (4)

and therefore similar to [5]. Here, we assume a statistical independence of all sensors and every single component refers to a probabilistic sensor model. The barometer information is evaluated using p(o_{t} | q_{t})_{baro}, whereby absolute position information is given by p(o_{t} | q_{t})_{ib} for iBeacons and by p(o_{t} | q_{t})_{wifi} for Wi-Fi.

A. Evaluation

The probability of currently residing on a floor is evaluated using the smartphone’s barometer. Environmental influences are circumvented by using relative pressure values instead of absolute ones. To reduce the impact of noisy sensors, we calculate the average \overline{ρ_{rel}} of several sensor readings, carried out while the pedestrian chooses his destination. This average serves as relative base for all future measurements. Likewise, we estimate the sensor’s uncertainty σ_{baro} for later use within the evaluation step.

In order to evaluate the relative pressure readings, we need a prediction to compare them with. Therefore, each transition from q_{t−1} to q_{t} estimates the state’s relative pressure prediction ρ̂_{rel} by tracking every height-change (z-axis):

qρ̂_{rel} ; t = qρ̂_{rel} ; t−1 + Δz · b, Δz = q_{t−1}^{z} − q_{t}^{z}, b ∈ ℝ . (5)

In eq. (5), b denotes the common pressure change in \frac{hPa}{m}. The evaluation step for time t compares every predicted relative pressure qρ̂_{rel} ; t with the observed one o_{t}^{ρ_{rel}} using a normal distribution with the previously estimated σ_{baro}:

p(o_{t} | q_{t})_{baro} = N (o_{t}^{ρ_{rel}} | qρ̂_{rel} ; t, σ_{baro}^{2}) .  (6)

The smartphone’s Wi-Fi and iBeacon component provides an absolute location estimation by measuring the signalstrengths of nearby transmitters. The positions of detected access-points (AP) and iBeacons are known beforehand. Using the wall-attenuation-factor signal strength prediction model [1], we are able to compare each measurement with a corresponding estimation. To infer this estimation, the prediction model uses the 3D distance d and the number of floors Δf between transmitter and the state-in-question q:

P_{r}(d, Δf) = P_{0} − 10γ log_{10} \frac{d}{d_{0}} + Δfβ ,  (7)

In eq. (7), there are three more parameters per AP. The signalstrength P_{0} measurable at a distance d_{0} (usually 1 m), a pathloss exponent γ describing the transmitter’s environment and the attenuation per floor β. To reduce the system’s setup time, we use the same three values for all access-points at the cost of accuracy. All parameters are chosen empirically. Further details on how to determine this parameters exactly, can be found in [18].

The same holds for the iBeacon component, except P_{0}, which is broadcasted by each beacon. However, as iBeacons cover only a small area, γ is usually much smaller compared to the one needed for Wi-Fi.

As transmitters are assumed to be statistically independent, the overall probability to measure their predictions at a given location is:

p(o_{t} | q_{t})_{wifi} = ∏_{i=1}^{n} N (s_{wifi}^{i} | P_{r}(d_{i}, Δf_{i}), σ_{wifi}^{2}) . (8)

B. Transition

The transition-distribution p(q_{t} | q_{t−1}) is sampled via random walks on a graph G = (V, E), which is generated from the buildings floorplan [5]. Each walk starts at q_{t−1} and uses adjacent edges e_{i,j} connecting two vertices v_{i}, v_{j} ∈ V until a certain distance d_{step} is reached. Hereby, the position of any q is represented by the position pos(v_{i}) of its corresponding vertex v_{i}. This approach draws only valid movements, as ambient conditions (walls, doors, stairs, etc.) are considered.

While sampling, to-be-walked edges are not chosen uniformly, but depending on a probability p(e_{i,j}). The latter depends on several constraints and recent sensor-readings from the smartphone. Using those readings directly within the transition step provides a more robust posterior distribution. Adding them to the evaluation instead, would lead to sample impoverishment due to the used Monte Carlo methods [19].

We assume the pedestrian’s desired destination to be known beforehand. This prior knowledge is incorporated during the random walk using p(e_{i,j})_{path}, which is a simple heuristic, favouring movements (edges) approaching his chosen destination with a ratio of 0.9 : 0.1 over those, departing from the destination [5]. The underlying shortest-path uses Dijkstra’s algorithm with special weight (distance) metric, considering special architectural facts:

Normally, the shortest-path calculated for a narrow grid would stick unnaturally close to obstacles like walls. To ensure realistic (human like) path estimations, we include architectural knowledge within Dijkstra’s edge-weight function [5]: Each vertex’s distance from the nearest wall is used to artificially increase the edge-weight and thus prevent the shortest-path from clinging to walls. Obviously this has a negative effect on doors which are surrounded by walls. Therefore, doors are detected and favoured by decreasing their edge-weight.

Steps and turns are detected using the smartphone’s IMU, implemented as described in [1]. The number of steps detected since the last transition is used to estimate the to-be-walked distance d_{step} by assuming a fixed step-size with some deviation:

d_{step} = o_{t−1}^{n_{steps}} · s_{step} + N (0, σ_{d_{step}}^{2}) .  (9)

Turn-Detection supplies the magnitude of the detected heading change by integrating over the gyroscope’s change since the last transition. Together with some deviation and the state’s previous heading, the magnitude is used to estimate the current state’s heading:

θ_{walk} = q_{t}^{θ} = q_{t−1}^{θ} + o_{t−1}^{Δθ} + N (0, σ_{θ_{walk}}^{2})  (10)

During the random walk, edges should satisfy the current heading and are thus drawn according to their resemblance:

p(e_{i,j})_{head} = p(e_{i,j} | θ_{walk}) = N (∠e_{i,j} | θ_{walk}, σ_{head}^{2}) . (11)

While the distribution eq. (11) does not integrate to 1.0 due to circularity of angular data, in our case, the normal distribution can be assumed as sufficient for small enough σ^{2}.

Additionally we perform a simple activity detection for the pedestrian, able to distinguish between several actions Ω ∈ {unknown, standing, walking, stairs up, stairs down}. For this, the sensor signals are split in sliding windows. Each window has a length of one second and overlaps 500 ms with its prior window. We use a naive Bayes classifier with two features. The first one is the variance of the accelerometer’s magnitude within a window. The second feature is the difference between the last and first barometer measurement of the particular window. Based on these features the classifier assigns an activity to each of the sliding windows. Similarly to the above, this knowledge is then evaluated when walking the grid: Edges e_{i,j} matching the currently detected activity are favoured using p(e_{i,j})_{act} = 0.8 and 0.2 otherwise. If no information of the current activitiy could be obtained, no influence is exerted on the edges.

IV. Smoothing

The main purpose of this work is to provide MC smoothing methods in context of indoor localisation. As mentioned before, those algorithms are able to compute probability distributions in the form of p(q_{t} | o_{1:T}) and are therefore able to make use of future observations between t and T, where t << T. In the following we discuss the algorithmic details of the forward-backward smoother and the backward simulation. Further, a novel approach for incorporating them into the localisation system is shown.

Algorithm 1 Forward-Backward Smoother 1: for t = 1 to T do ▷ Filtering 2: Obtain the weighted trajectories {W_{t}^{i}, X_{t}^{i}}_{i=1}^{N}

TODO Filtering hier genauer beschreiben?

3: end for 4: for i = 1 to N do ▷ Initialization 5: Set W_{T |T}^{i} = W_{T}^{i} 6: end for 7: for t = T − 1 to 1 do ▷ Smoothing 8: for i = 1 to N do

9: W_{t|T}^{i} = W_{t}^{i} [∑_{j=1}^{N} W_{t+1|T}^{j} \frac{p(X_{t+1}^{j}|X_{t}^{i})}{∑_{k=1}^{N} W_{t}^{k} p(X_{t+1}^{j}|X_{t}^{k})}]

10: end for 11: end for

A. Forward-backward Smoother

The forward-backward smoother (FBS) of [20] is a well established alternative to the simple filter-smoother. The foundation of this algorithm was again laid by Kitagawa in [21]. An approximation is given by

p(q_{t} | o_{1:T}) ≈ ∑_{i=1}^{N} W_{t|T}^{i} δ_{X}_{t}^{i}(q_{t}) ,  (12)

where p(q_{t} | o_{1:T}) has the same support as the filtering distribution p(q_{t} | o_{1:t}), but the weights are different. This means, that the FBS maintains the original particle locations and just reweights the particles to obtain a smoothed density. δ_{X}_{t}^{i} denotes the Dirac delta function. The complete FBS can be seen in algorithm IV-A in pseudo-algorithmic form. At first, the algorithm obtains the filtered distribution (particles) by deploying a forward step at each time t. Then the backward step for determining the smoothing distribution is carried out. The weights are obtained through the backward recursion in line 9.

By reweighting the filter particles, the FBS improves the simple filter-smoother by removing its dependence on the inheritance (smoothed) paths [22]. However, by looking at algorithm IV-A it can easily be seen that this approach computes in O(N^{2}), where the calculation of each particle’s weight is an O(N) operation. To reduce this computational bottleneck, [23] introduced a solution using algorithms from N-body simulation. By integrating dual tree recursions and fast multipole techniques with the FBS, a run-time cost of O(N log N) can be achieved.

B. Backward Simulation

For smoothing applications with a high number of particles, it is often not necessary to use all particles for smoothing. This decision can for example be made due to a high sample impoverishment and/or highly accurate sensors. By choosing a good sub-set for representing the posterior distribution, it is theoretically possible to further improve the estimation.

Therefore, [10] presented the backward simulation (BS). Where a number of independent sample realisations from Algorithm 2 Backward Simulation Smoothing 1: for t = 1 to T do ▷ Filtering 2: Obtain the weighted trajectories {W_{t}^{i}, X_{t}^{i}}_{i=1}^{N} 3: end for 4: for k = 1 to N_{sample} do 5: Choose q̃_{T}^{k} = X_{T}^{i} with probability W_{T}^{i} ▷ Initialize 6: for t = T − 1 to 1 do ▷ Smoothing 7: for j = 1 to N do 8: W_{t|t+1}^{j} = W_{t}^{j} p(q̃_{t+1} | X_{t}^{j}) 9: end for 10: Choose q̃_{t}^{k} = X_{t}^{j} with probability W_{t|t+1}^{j} 11: end for 12: q̃_{1:T}^{k} = (q̃_{1}^{k}, q̃_{2}^{k}, ..., q̃_{T}^{k}) is one approximate realisation

from p(q_{1:T} | o_{1:T}) 13: end for

the entire smoothing density are used to approximate the smoothing distribution. This method can be seen in algorithm 2 in pseudo-algorithmic form. Again, a particle filter is performed at first and then the smoothing procedure gets applied. Here, q̃_{t} is a random sample drawn approximately from p(q_{t} | q̃_{t+1}, o_{1:T}). For example q̃_{t} could be chosen by selecting particles within a cumulative frequency. Therefore q̃_{1:T} = (q̃_{1}, q̃_{2}, ..., q̃_{T}) is one particular sample realisation from p(q_{1:T} | o_{1:T}). Further independent realisations are obtained by repeating the algorithm until the desired number N_{sample} is reached. The computational complexity for one particular realisation is O(N). However, the computations are then repeated for each realisation drawn [10].

C. Transition for Smoothing

As seen above, both algorithms are reweighting particles based on a state transition model. Unlike the transition presented in section III-B, it is not possible to just draw a set of new samples. Here, p(q_{t+1} | q_{t}) needs to provide the probability of the known future state q_{t+1} under the condition of the current state q_{t}. In case of indoor localisation using particle filtering, it is necessary to not only provide the probability of moving to a particle’s position under the condition of its ancestor, but also of all other particles at time t. The smoothing transition model therefore calculates the probability of being in a state q_{t+1} in regard to previous states and the pedestrian’s walking behaviour. This means that a state q_{t} is more likely if it is a proper ancestor (realistic previous position) of a future state q_{t+1}. In the following a simple and inexpensive approach for receiving this information will be described.

By writing

p(q_{t+1} | q_{t}, o_{t})_{step} = N (Δd_{t} | μ_{step}, σ_{step}^{2})  (13)

we receive a statement about how likely it is to cover a distance Δd_{t} between two states q_{t+1} and q_{t}. In the easiest case, Δd_{t} is the linear distance between two states. Of course, based on the graph structure, one could calculate the shortest path between both and sum up the respective edge lengths. However, this requires tremendous calculation time for negligible improvements. Therefore this is not further discussed within this work. The average step length μ_{step} is based on the pedestrian’s walking speed and σ_{step}^{2} denotes the step length’s variance. Both values are chosen depending on the activity Ω recognized at time t. For example μ_{step} gets smaller while a pedestrian is walking upstairs, than just walking straight. This requires to extend the smoothing transition by the current observation o_{t}. Since q is hidden and the Markov property is satisfied, we are able to do so.

The heading information is incorporated using

p(q_{t+1} | q_{t}, o_{t})_{turn} = N (Δα_{t} | Δθ, σ_{turn}^{2}) ,  (14)

where Δα_{t} is the absolute angle between q_{t+1} and q_{t} in the range of [0, π]. The relative angular change Δθ is then used to receive a statement about how likely it is to walk in that particular direction. Again the normal distribution of eq. (14) does not integrate to 1.0. Therefore the same assumption as in eq. (11) has to be made.

To further improve the results, especially in 3D environments, the vertical (non-absolute) distance Δz between two successive states is used as follows:

p(q_{t+1} | q_{t}, o_{t})_{baro} = N (Δz | μ_{z}, σ_{z}^{2}) .  (15)

This assigns a low probability to false detected or misguided floor changes. Similar to eq. (13) we set μ_{z} and σ_{z}^{2} based on the activity recognised at time t. Therefore, μ_{z} is the expected change in z-direction between two time steps. This means, if the pedestrian is walking alongside a corridor, we set μ_{z} = 0. In contrast, μ_{z} is positive while walking downstairs or otherwise negative for moving upstairs. The size of μ_{z} and also μ_{step} could be a predefined value or set dynamically based on the measured vertical and linear acceleration.

Looking at eq. (13) to eq. (15), obvious similarities to a sensor fusion process can be seen. By assuming statistical independence between those three, the probability density of the smoothing transition is given by

p(q_{t+1} | q_{t}) = p(q_{t+1} | q_{t}, o_{t})_{step} p(q_{t+1} | q_{t}, o_{t})_{turn} p(q_{t+1} | q_{t}, o_{t})_{baro} .  (16)

It is important to notice, that all particles at each time step t of the forward filtering need to be saved. Therefore, the memory requirement increases proportional to the processing time.

V. Experiments

The experiments were carried out on all floors (0 to 3) of the faculty building. Each floor is about 77 m x 55 m in size, with a ceiling height of 3 m. To resemble real-world conditions, the evaluation took place during an in-house exhibition. Thus, many places were crowded and Wi-Fi signals attenuated. As can be seen in fig. 2 we arranged 4 distinct walks, covering different distances, critical sections and uncertain decisions leading to multimodalities. The ground truth is measured by recording a timestamp at marked spots on the walking route.

path 1 path 2 path 3 path 4 path 1 path 2 path 3 path 4

Fig. 2. The four paths that were part of the evaluation. Starting positions are marked with black circles. For a better visualisation they were slightly shifted to avoid overlapping.

When passing a marker, the pedestrian clicked a button on the smartphone application. Between two consecutive points, a constant movement speed is assumed. Thus, the ground truth might not be 100 % accurate, but fair enough for error measurements. The approximation error is then calculated by comparing the interpolated ground truth position with the current estimation. Especially in the context of smoothing, it is also very interesting to exclude the temporal delay from the error calculations and measure only the positional difference between estimated and ground truth path. This gives a statement about the extent to which the smoothed path superficially improves compared to the filtered one.

All walks start with a uniform distribution (random position and heading) as prior for q_{0}. To allow the system to stabilize its initial state, the first few estimations are omitted from error calculations. Even though, the error during the following few seconds is expected to be much higher than the error when starting with a well known initial position and heading.

The measurements were recorded using a Motorola Nexus 6 and a Samsung Galaxy S5. As the Galaxy’s Wi-Fi can not be limited to the 2.4 GHz band only, its scans take much longer than those of the Nexus: 3500 ms vs. 600 ms. Additionally, the Galaxy’s barometer sensor provides fare more inaccurate and less frequent readings than the Nexus does. This results in a better localisation using the Nexus smartphone. The computation for both filtering and smoothing was done offline using the aforementioned CONDENSATION algorithm and multinomal (cumulative) resampling. Unless explicitly stated, the state was estimated using the weighted arithmetic mean of the particles.

As mentioned earlier, the position of all access-points (about 5 per floor) is known beforehand. Due to legal terms, we are not allowed to depict their positions and therefore omit this information within the figures. Additionally, we used five iBeacons for slight enhancements in some areas. The empirically chosen values for Wi-Fi were P_{0_{wifi}} = −46 dBm, γ_{wifi} = 2.7 , β_{wifi} = 8 dB, and γ_{ib} = 1.5 for the iBeacons, respectively. The system was tested by omitting any time-consuming calibration processes for those values. We therefore expect the

unkown standing walking stairs up stairs down 0 s 20 s 40 s 60 s 80 s 100 s 120 s 140 s 1  2  3 4 5  6  7 8 9 10

Fig. 3. activity recognition

filtered  smoothed  0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007

Fig. 4. particles. green = avg50, black = avg. gnuplot zickt bei der legende

localisation process to perform generally worse compared to standard fingerprinting methods [24]. However, incorporating prior knowledge and smoothing will often compensate for those poorly chosen system parameters.

For the filtering we used σ_{wifi} = σ_{ib} = 8.0 as uncertainties, both growing with each measurement’s age. While the pressure change was assumed to be 0.105 \frac{hPa}{m}, all other barometerparameters are determined automatically (see III-A). The step size s_{step} for the transition was configured to be 70 cm with an allowed derivation of 10 % . The heading deviation was set to 25^{∘}. Edges departing from the pedestrian’s destination are downvoted using κ = 0.9. For smoothing we set σ_{turn}^{2} = 5^{∘} and σ_{z}^{2} = 0.25 cm. If walking or unknown are the current activities, μ_{step} = 0.7 m, μ_{step} = 0.5 and μ_{z} = 0.0 m are used. Walking upstairs sets μ_{step} = 0.4 m, σ_{step}^{2} = 0.2 and μ_{z} = −0.3 m, otherwise μ_{step} = 0.5 m, σ_{step}^{2} = 0.3 and μ_{z} = 0.3 m for walking downstairs.

By adding the activity recognition to the system of [5], we are able to further improve the overall localisation results. The approximation error decreases by an average of XX cm for all 4 paths on 10 MC runs. Due to this additional knowledge, the state transition samples mostly depending upon the current activity and therefore limits the possibility of false floor changes. Fig. 3 shows the recognized activities for path 4 using the Nexus 6. Despite a short misdetection caused by faulty pressure readings, the recognition can be considered to be very robust and accurate.

At first, both FBS and BS are compared in context of fixedinterval smoothing. As a reminder, fixed-interval smoother are using all observations until time T and therefore run offline, after the filtering procedure is finished. Thus, we calculate only the positional error between estimation and ground truth, since timely information are negligible. In contrast to BS, the FBS

smoothed (FBS) smoothed (BS) smoothed (FBS) smoothed (BS) 2 3  4 5  6  7 89 10

Fig. 5. Comparison between FBS (red) and BS (blue) on path 4 (black). Both were approximated using 2500 particles and 500 sample realisations for BS. The measurements were recorded using the Nexus 6. TODO Filter mit rein

filtered smoothed (BS) filtered smoothed (BS)

Fig. 6. A situation where BS smoothing (blue) was not able to improve the filtering results (green). Two main factors are causing this: an initial position within a detached room and inaccurate pressure readings given by the Galaxy S5.

is not able to improve the results using the weighted arithmetic mean for estimating the current position. As illustrated in fig. 4, the estimated position (black dot) for filtering and FBS is identical although the weights are highly different. To address this problem, we are calculating the average state using the 50 % best weighted particles (green dot). For evaluating the FBS this estimation method was applied to filtering and smoothing likewise. We deployed 10 MC runs using 2500 particles for approximation. Now, the positional average error along all 4 paths could be improved from 2.49 m to 2.2 m for the Galaxy and from 1.76 m to 1.54 m for the Nexus.

However, BS still outperforms the FBS with an average error of 1.74 m for the Galaxy and 1.41 m for the Nexus on all 4 paths using the same number of particles and 500 sample realisations. A visual example comparing both smoothing methods on path 4 is illustrated in fig. 5. The estimation of BS looks way more realistic and adapts better to the ground truth path. However, in this particular example the FBS starts at an earlier position (cf. fig. 5 seg. 2), better handling the initial uniform distribution. Another advantage of BS over FBS, is the ability to still improve the filtering results even while reducing the number of particles radical. For example 50 particles and 25 sample realisations are providing reliable estimations similar to above experiments, though the risk of losing track is higher.

Despite the very good outcomes provided by both interval smoother, there are some rare situations in which smoothing does not improve the filtered estimation or even improves

0.0 m 2.0 m 4.0 m 6.0 m 8.0 m0 s  20 s  40 s  60 s  80 s  100 s  120 s 2 3  4  5  6  7  8  9 10 filtered  BS  FBS

Fig. 7. Hier muss ich noch was schreiben. Bin aber seeeehr muede. Geh jetzt ins Bett. Nachti.

the visual path. For example fig. 6 depicts such a situation for path 3 using BS and measurements provided by the Galaxy S5. Here, the estimation was not able to change floors correctly due to faulty pressure readings. Additionally, the initial position was located within a detached room. This shows that the smoothing results are of course highly depend upon the filtering performance.

At next, we discuss the advantages and disadvantages of utilizing FBS and BS as fixed-lag smoother. Compared to fixed-interval smoothing, timely errors are now of higher importance due to an interest on real-time localization. Especially interesting in this context are small lags τ < 10 considering filter updates near 500 ms. Fig. 7 illustrates the different approximation errors alongside path 4 using 500 particles, 100 sample realisations for BS and a fixed-lag τ = 5. Due to the small number of sample realisations for BS and the different estimation method for FBS, the errors are changing very frequently in contrast to the filter. It can be seen that again BS provides a better overall estimation, especially in areas where the user is changing floors (cf. fig. 7 seg. 4, 9). Besides the positional quality, also the timely error could be reduced by both algorithms in fig. 7. Once more, the BS outperforms the FBS by providing an overall approximation error of 4.86 m for the Galaxy and 2.97 m for the Nexus by filtering with 5.68 m and 3.15 m respectively. Whereas FBS is not able to reduce the filtering error, obtaining 5.59 m using the Nexus and 3.12 m the Galaxy. However, this does not mean, that the FBS is unpracticable for problems of fixed-lag smoothing. By looking at the data in detail, high errors similar as seen in fig. 7 occur on path 3, which cause the estimation to behave more natural instead of walking the supposed path. This phenomena could be observed for both smoothers respectively. Also note the difference between this error, including timely information, and the positional error used before. Similar to fixed-interval smoothing, decreasing the number of particles does not necessarily worsen the estimation. In most cases smoothing compensates for this reduction and maintains the good results.

Besides changing the number of particles, it is also possible the variate the lag. As one would expect, increasing the lag causes the smoothed estimation to approach the results provided by fixed-interval smoothing. This can be verified by looking at fig. ??, which is a detailed view of segment XX of path 4 (cf. fig. 5). It is obvious that a lag of 30 time steps has access to much more future observations and is therefore able to obtain such a result. Considering an update interval of 500 ms, a lag of 30 would however mean that the smoother is 15 s behind the filter. Nevertheless, there are practical applications like accurately verifying hit checkpoints or continuously optimizing a recurring segment of the path.

VI. Conclusion Within this work a novel approach for utilizing the forwardbackward smoother and backward simulation to problems of indoor localisation was presented. Both were implemented as fixed-lag and fixed-interval smoother. It was shown that smoothing methods are able to decrease the estimation error and improving the overall localisation. Especially fixed-lag smoothing is a great tool for runtime support by reducing timely errors and improving the overall estimation with affordable costs. However, a fixed-lag smoother is not able to change the lag dynamically, as its name suggests. Therefore, a dynamic-lag smoother could be able to further improve the estimation by considering higher lags in critical areas.

Finally, the smoothing transition does not use any information provided by the underlying graph structure. This would allow to use environmental informations and to replace the current line-of-sight model with a graph-based one. By incorporating the Wi-Fi’s signal strength measurements a more advanced smoothing transition should be able to compensate for faulty Wi-Fi measurements and the hereby resulting jumps between positions.

References

[1] F. Ebner, T. Fetzer, L. Köping, M. Grzegorzek, and F. Deinzer, “Multi Sensor 3D Indoor Localisation,” in Indoor Positioning and Indoor Navigation (IPIN), Int. Conf. on, 2015, pp. 1–10.

[2] H. Nurminen, M. Koivisto, S. Ali-Loytty, and R. Piche, “Motion Model for Positioning with Graph-based Indoor Map,” in Int. Conf. on Indoor Positioning and Indoor Navigation (IPIN). IEEE, 2014, pp. 646–655.

[3] S. Hilsenbeck, D. Bobkov, G. Schroth, R. Huitl, and E. Steinbach, “Graph-based Data Fusion of Pedometer and WiFi Measurements for Mobile Indoor Positioning,” in Proc. of the 2014 ACM Int. Joint Conf. on Pervasive and Ubiquitous Computing - UbiComp ’14 Adjunct. New York, NY, USA: ACM Press, 2014, pp. 147–158.

[4] B. Li, B. Harvey, and T. Gallagher, “Using barometers to determine the height for indoor positioning,” in Indoor Positioning and Indoor Navigation (IPIN), Int. Conf. on, 2013, pp. 1–7.

[5] F. Ebner, T. Fetzer, M. Grzegorzek, and F. Deinzer, “On Prior Navigation Knowledge in Multi Sensor Indoor Localisation,” in On Review, 2016, pp. 1–8.

[6] S. Thrun, W. Burgard, and D. Fox, Probabilistic Robotics. Cambridge, MA, USA: The MIT Press, 2005.

[7] A. Doucet and A. M. Johansen, “A Tutorial on Particle Filtering and Smoothing: Fifteen years later,” in The Oxford Handbook of Nonlinear Filtering, D. Crisan and B. Rozovsky, Eds. Oxford University Press, 2011.

[8] Z. Chen, “Bayesian filtering: From kalman filters to particle filters, and beyond,” Statistics, vol. 182, no. 1, pp. 1–69, 2003.

[9] A. Doucet, S. Godsill, and C. Andrieu, “On sequential monte carlo sampling methods for bayesian filtering,” Statistics and Computing, vol. 10, no. 3, pp. 197–208, 2000.

[10] S. J. Godsill, A. Doucet, and M. West, “Monte carlo smoothing for nonlinear time series,” Journal of the American Statistical Association, vol. 99, no. 465, pp. 156–168, 2004.

[11] G. Kitagawa, “Monte carlo filter and smoother for non-gaussian nonlinear state space models,” Journal of Computational and Graphical Statistics, vol. 5, no. 1, pp. 1–25, 1996.

[12] M. Isard and A. Blake, “A Smoothing Filter for Condensation,” in Computer VisionECCV’98. Springer, 1998, pp. 767–781.

[13] P. Pérez, J. Vermaak, and A. Blake, “Data Fusion for Visual Tracking with Particles,” Proceedings of the IEEE, vol. 92, no. 3, pp. 495–513, 2004.

[14] E.-S. Platzer, F. Deinzer, D. Paulus, and J. Denzler, “3D Blood Flow Reconstruction from 2D Angiograms,” in Bildverarbeitung für die Medizin 2008, ser. Informatik Aktuell, T. Tolxdorff, J. Braun, T. M. Deserno, A. Horsch, H. Handels, and H.-P. Meinzer, Eds. Springer Berlin Heidelberg, 2008, pp. 288–292.

[15] Y. Hu, P. Baraldi, F. Di Maio, and E. Zio, “A Particle filtering and Kernel Smoothing-Based Approach for New Design Component Prognostics,” Reliability Engineering and System Safety, vol. 134, pp. 19–31, 2014.

[16] A. S. Paul and E. A. Wan, “RSSI-Based Indoor Localization and Tracking using Sigma-Point Kalman Smoothers,” IEEE Journal on Selected Topics in Signal Processing, vol. 3, no. 5, pp. 860–873, 2009.

[17] L. Köping, M. Grzegorzek, and F. Deinzer, “Probabilistic Step and Turn Detection in Indoor Localization,” in Conf. on Data Fusion and Target Tracking 2014: Algorithms and Applications (DFTT 2014), Liverpool, UK, 2014, pp. 1–7.

[18] S. Y. Seidel and T. S. Rappaport, “914 MHz path loss prediction models for indoor wireless communications in multifloored buildings,” IEEE Transactions on Antennas and Propagation, vol. 40, no. 2, pp. 207– 217, 1992.

[19] M. Isard and A. Blake, “CONDENSATION - Conditional Density Propagation for Visual Tracking,” Int. Journal of Computer Vision, vol. 29, no. 1, pp. 5–28, 1998.

[20] A. Doucet, S. Godsill, and C. Andrieu, “On Sequential Monte Carlo Sampling Methods for Bayesian Filtering,” Statistics and Computing, vol. 10, no. 3, pp. 197–208, 2000.

[21] G. Kitagawa, “Non-gaussian state-space modeling of nonstationary time series,” Journal of the American Statistical Association, vol. 82, no. 400, pp. 1032–1041, 1987.

[22] P. Fearnhead, D. Wyncoll, and J. Tawn, “A sequential smoothing algorithm with linear computational cost,” Biometrika, vol. 97, no. 2, pp. 447–464, 2010.

[23] M. Klaas, M. Briers, N. De Freitas, A. Doucet, S. Maskell, and D. Lang, “Fast particle smoothing: If i had a million particles,” in Proc. of the 23rd Int. Conf. on Machine learning. ACM, 2006, pp. 481–488.

[24] V. Honkavirta, T. Perl, S. Ali-Lytty, and R. Piche, “A comparative survey of wlan location fingerprinting methods,” Workshop on Positioning, Navigation and Communication, vol. 6, no. 1, pp. 243–251, 2009.