Simultaneous localization and mapping in wireless sensor networks

Mobile device localization in wireless sensor networks is a challenging task. It has already been addressed when the WiFi propagation maps of the access points are modeled deterministically or estimated using an offline human training calibration. However, these techniques do not take into account the environmental dynamics. In this paper, the maps are assumed to be made of an average indoor propagation model combined with a perturbation field which represents the influence of the environment. This perturbation field is embedded with a distribution describing the prior knowledge about the environmental influence. The device is localized with Sequential Monte Carlo methods and relies on the estimation of the propagation maps. This inference task is performed online, using the observations sequentially, with a new online Expectation Maximization based algorithm. The performance of the algorithm is illustrated with Monte Carlo experiments using both simulated data and a true data set. & 2014


Introduction
In this paper, we consider a WiFi communication network made up of a mobile device, a server and WiFi access points (AP).In this context, a key step to localize the mobile device is to estimate the WiFi signal strength at different positions in the environment.However, in an indoor environment, signals may experience complex attenuation such as shadowing or reflection.
Different techniques can be used to approximate the WiFi propagation map of each AP.In Gorce et al. [19], a deterministic model based on the characteristics of the surrounding AP and on the localization of the obstacles in the environment is introduced.In Bahl and Padmanabhan [2] and Evennou and Marx [16], a previous offline training phase is performed.In this site survey, the signal strength indicator (RSSI) received from different AP is measured at some previously determined positions.This allows us to build an accurate estimation of the signal strength, but only for a finite number of points.Ferris et al. [17] provide a method to extend these measures to the entire map using Gaussian processes techniques.
In this paper, we propose an estimation method that does not require any calibration procedure.The propagation maps are estimated online, without storing the observations, using the data sent by the mobile device.Any modification in the signal propagation (due to new obstacles for instance) affects the data sent by the mobile device.Then, while these changes deteriorate the accuracy of localization systems based on fixed estimators of the propagation maps, our estimation procedure takes them into account.Thus, as illustrated in Section 5.2, the accuracy of our localization method improves with time instead of degrading.
We use a semiparametric statistical model introduced in [15, Chapter 5]: the propagation maps are made of a parametric average indoor model and a nonparametric perturbation field.This model combines a prior knowledge on signal propagation with random perturbations due to the obstacles.Based on the data collected by the mobile device, the parameters and the perturbation field can be estimated.The proposed procedure relies on a new online Expectation Maximization (EM) based algorithm introduced in Le Corff and Fort [21,22] and on Sequential Monte Carlo methods.The device position can be simultaneously estimated using the weighted samples produced by the Sequential Monte Carlo method.
The structure of this paper is the following.Section 2 details the approach of the paper in comparison with other methods for localization and mapping in wireless sensor networks.Section 3 describes the model and defines the notations.Section 4 introduces our algorithm for the considered Simultaneous Localization and Mapping (SLAM) problem.Section 5 illustrates the performance of this algorithm with numerical experiments.

Wireless sensor networks
In a radio network field, the signal strength is measured by the RSS (received signal strength, in dBm).Each WiFi connected device can compute the RSS as it is needed to associate the device with the AP providing the best signal/ noise ratio.Localization systems based on RSS measurements allow us to locate any WiFi connected device.
Remark 1. Different devices might have different methods to compute the RSS.It is common to use the terminology of RSSI (received signal strength indicator) expressed without unit to name the information on the signal level provided by a WiFi device.To overcome the issue of information disparity between WiFi connected devices, a previous RSSI to RSS conversion might be needed for each localized device.The conversion rule is specific to the device's WiFi card.For the sake of clarity, we assume in this paper that the mobile device's RSSI corresponds to the standard RSS.
Using RSS to estimate the position of the device is challenging.As represented in Fig. 1, RSS is highly unstable as it fluctuates considerably between two consecutive measures at the same position.It is common to describe the RSS variations using a Gaussian representation.
Despite the instant variations of the RSS, its mean value strongly depends on the position of the device.The function that returns the expected RSS at each position is called the propagation map.In free spaces, signals propagate in straight line from the emitter to the receptor.The propagation maps are isotropic and can be described using few parameters, see Friis [18].Therefore, there exist two parameters c and d such that the strength of a signal received at x and emitted at O is c þ d logð J x À OJ Þ where c and d depend on the network characteristics.
On the other hand, indoor propagation of WiFi signals is not isotropic.Fig. 2 represents the propagation map of a WiFi signal in an indoor environment.This figure was built deterministically by Gorce et al. [19] using the physical properties of electromagnetic signals.

General model for signal propagation
Let X k be the device position at time k.The received signal strength vector measured by the device, denoted by Y k , is written as In the following, the superscript ⋆ is used to name the true value of every unknown parameter involved in the model.fϵ k g k Z 0 are i.i.d.multidimensional random variables with distribution N ð0; s ⋆;2 IdÞ where Id is the identity matrix.Each component of F ⋆ k ðxÞ represents the expected RSS at position x relative to an AP.The time dependency of F ⋆ brings the environmental effects on the propagation maps into relief.The distribution of fϵ k g k Z 0 implies that the noise affecting the RSS of different AP are independent.This assumption is somehow strong, however, the correlation between the RSS of different AP can be hardly taken into account and can strongly depend on the environment configuration.

Contribution of the paper and comparison with the state of the art
Indoor localization requires to overcome two challenges.One has to design a good approximation of F ⋆ k and then to use this approximation to estimate the position X k corresponding to a given measure Y k .
The contribution of this paper is a new estimation procedure of the propagation maps.This new method is combined to a localization algorithm to highlight the relevance of our estimation procedure.Using our propagation maps estimator, many positioning algorithms may be considered.We do not address a comparison of the accuracy between our method and the state of art but a new way to improve indoor localization deployment at a large scale.
F ⋆ k can be approximated deterministically using physical properties (see [19,26]).Such constructions need a precise description of the environment such as the position and the composition of the walls and furniture present in the environment.They also need a fine knowledge about the different factors influencing signal propagation such as the humidity level in the environment.Without this precise description, the propagation maps built using deterministic methods may be far from reality.
Most of the existing WiFi based localization systems use a preliminary measurement campaign (offline) which is also called fingerprinting technique.A human operator performs a site survey by measuring the RSS from different AP at some fixed positions in the environment.This set of measures, associated with positions, can be directly stored in a database [2] to be used in the localization system as a reference.These measures can also be extrapolated to the entire map using kriging methods such as in Ferris et al. [17].Users must compile a fairly dense radio map comprising many RSS measurements at many sampled points to attain reasonable positioning accuracy.Moreover, these methods suppose that the propagation maps remain constant which, with our notation, means that the propagation function F ⋆ does not depend on time.The measurement campaign can be seen as an instant "photograph" of the propagation map and has to be regularly performed in order to maintain the accuracy of the system.These problems have particularly been spotted by Chen et al. [9] and Madigan et al. [24].Chen et al. [9] introduce RFID sensors in the environment to perform passive site survey.Madigan et al. [24] use a hierarchical Bayesian model to localize the mobile without site survey but rely on isotropic propagation maps which might lead to a bad estimation in complex indoor environments (although they study a more elaborate model that includes "corridors effects").
In this paper, we present a new estimation procedure of the propagation maps which does not involve any measurement campaign or additional sensors.The considered model does not assume any knowledge on the environment apart from the position of the AP.The propagation maps are estimated using an online Expectation Maximization (EM) based algorithm using the RSS measurements collected by the device.A similar approach can be found in Chai and Yang [8] which uses the classical EM algorithm to refine the propagation maps estimators obtained using a preliminary site survey.The first substantial benefit of our method concerns its ability to be implemented in any building.Moreover, unlike fingerprinting based methods, the precision of our localization system does not degrade with time as each measure Y k is used to improve the propagation maps estimators.The computation of the estimators requires sufficient statistics.These statistics summarize the information contained in all the past measures since they were last reset.If there are regular environmental modifications that strongly affect the WiFi propagation, regularly resetting these sufficient statistics is enough to improve the estimation.Then, our system is more robust than site survey based methods using static propagation maps estimators.
Once the expected RSS has been estimated for the whole map or for some fixed positions, the device can be localized.Bahl and Padmanabhan [2] use the nearest neighbor algorithm.With this algorithm, the strong variability of the RSS leads to a very unstable localization.We use particle filtering to track the mobile device.Such filters have already been used by Ferris et al. [17] in a similar way.Evennou and Marx [16] introduce particle filtering on the Voronoi diagram.Such filters provide a much more stable sequence of positions despite the high variability of the data.
Remark 2. There is no chance to identify fF ⋆ k g k Z 0 with the observations fY k g k Z 0 only.In the next section, we omit the time dependency of F ⋆ which might seem contradictory with our introduction.However, as stated in the above section, regularly resetting the sufficient statistics allows us to adapt the algorithm to environmental changes.

Model and assumptions
Let X k be the cartesian coordinates of the mobile device at time k in a two-dimensional compact space.This compact space represents the map of the one-floor building where the localization is performed.This continuous environment is discretized into a finite grid map, denoted by C. It is assumed that fX k g k Z 1 is a Markov chain taking values in C with initial distribution ν and Markov transition matrix given, for all ðx;  where J Á J denotes the usual euclidean norm in R 2 (the associated inner product is denoted by 〈Á; Á〉). a A R ⋆ þ depends on the average speed of the mobile and is assumed to be known.Therefore, the initial state X 0 is distributed according to ν and, for any k Z 1 and any Let B be the number of AP and jCj be the cardinality of C. In the sequel, for any B Â R jCj matrix A, we use the shorthand notation A j for the vector fA j;x g x A C and A j 2 for the vector fA 2 j;x g x A C .At each time step k Z1 and for each j A f1; …; Bg, the mobile device measures and sends to the server the observation Y k;j given by x is the j-th average indoor propagation term at position x.For all x A C and all jA f1; …; Bg, μ ⋆ j;x only depends on the distance between x and O j where O j is the known position of the j-th AP.In the sequel, we use the so-called Friis transmission equation given by Friis [18]: where c ⋆ j and d ⋆ j are parameters depending on the environment and where log is the logarithm to the base e. δ ⋆ j is an additive term due to random perturbations.A similar model of WiFi propagation maps using Gaussian processes can be found in Ferris et al. [17].It is assumed that the parameters fδ ⋆ j g B j ¼ 1 are embedded with the prior distribution π given, for any δ A R BjCj , by ; where for any matrix A, A T is the transpose of A and where Σ is assumed to be known.
, with mean 0 and covariance matrix s ⋆;2 I B , where I B is the identity matrix of size B Â B.
j will be referred to as the true propagation map of the j-th AP and as the true propagation maps.Fig. 3 displays a realization of δ ⋆ j (sampled from π) and the functions μ ⋆ j and F ⋆ j , defined on the grid C ¼ f0; …; 30g Â f0; …; 30g.The parameters used in this figure are O j ¼ ð15; 15Þ, and c ⋆ j , d ⋆ j and Σ are given in Section 5.In the sequel, we write where For any x A C and any k Z 1, the distribution of Y k conditionally to X k ¼ x has a density with respect to the Lebesgue measure on R B given by In the sequel, we aim at simultaneously estimating the device position and θ ⋆ using the observations fY k g k Z 1 .For any positive integer n, any observation set ðy 1 ; …; y n Þ shortly denoted by y 1:n , and any parameter θ ¼ ðc; d; δ; s 2 Þ, the likelihood of the observations y 1:n is given by Let n be a positive integer and Y 1:n be a set of observations.The estimator of θ ⋆ is set as one maximizer of the function: 4. Online estimation procedure

EM based algorithms to estimate the propagation maps
The EM algorithm is a well-known iterative algorithm to perform maximum likelihood estimation in hidden Markov models [13].An EM based algorithm can be introduced to maximize (5).Each iteration of this algorithm is decomposed into an E-step where the expectation of the complete data log-likelihood (log of the joint distribution of the states and the observations) conditionally on the observations is computed; and a M-step which updates the parameter estimate.Let Y 1:n be a fixed set of observations and θ i be the current map estimate.
(i) The E-step amounts to computing the conditional expectation where log p θ ðX 1:n ; Y 1:n Þ is the complete data loglikelihood and E θ i ½ÁjY 1:n is the conditional expectation given Y 1:n when the map is θ i .(ii) The M-step computes the new value θ i þ 1 by choosing one of the maps θ maximizing θ ¼ ðc; d; δ; s 2 Þ↦Q θ i ðY 1:n ; θÞþn À 1 log πðδÞ: Define, for any ðx; yÞ A C Â R B and any j A f1; …; Bg, the vectors where 1 x 0 ðxÞ equals 1 if x ¼ x 0 and 0 otherwise.The constant a being known, the intermediate quantity associated with the model presented in Section 3 can be written, up to an additive constant, as where (the dependence on θ i , n and Y 1:n is dropped from the notation for better clarity) Therefore, by (7), it is enough to compute S 1 , S 2;j and S 3;j , 1 rj rB, to maximize the function θ ¼ ðc; d; δ; s 2 Þ↦Q θ i ðY 1:n ; θÞþn À 1 log πðδÞ.The detailed computations to solve this optimization problem and to obtain the new map estimate θ i þ 1 are given in Algorithm 1 (where 1 is the vector of size jCj where each entry equals 1 and, for any vector v, diagðvÞ is the diagonal matrix with diagonal given by v).

21:
end for 22: This two step process can be repeated until the likelihood does not improve significantly.However, when the observations are obtained sequentially, this algorithm does not produce a new estimate as new observations are received.The mobile device localization requires an online method which does not store the data and which frequently updates the propagation maps.
Online variants of the EM algorithm have been proposed to obtain map estimates each time a new observation is available.In the case of i.i.d.observations, Cappé and Moulines [6] proposed the first EM based online algorithm.This algorithm replaces the exact computation of the sufficient statistics S 1 , S 2;j and S 3;j by a stochastic approximation step.In the case of hidden Markov models, when both the Fig. 3. Example of δ ⋆ j (sampled from π) and of the functions μ ⋆ j and observations and the states take a finite number of values (resp.when the state-space is finite) an online EM-based algorithm was proposed by Mongillo and Denève [25] (resp.by [5]).These algorithms combine an online approximation of the filtering distributions of the hidden states and a stochastic approximation step to compute an online approximation of the sufficient statistics.This has been extended to the case of general state-space models with Sequential Monte Carlo algorithms in Cappé [4], Del Moral et al. [12] and Le Corff et al. [23].More recently, Le Corff and Fort [21,22] proposed the Block Online EM (BOEM) algorithm in which the estimate is kept fixed as block of observations is received and is updated at the end of each block.See also Andrieu et al. [1] for an overview of online estimation procedures using Sequential Monte Carlo methods.

Proposed algorithm for online localization in wireless sensor networks
This paper introduces an EM algorithm for online localization in wireless sensor networks based on the algorithm introduced in Le Corff and Fort [21,22].Let fτ k g k Z 1 be a sequence of block-sizes and define T 0 ¼ def 0 and, for any k Z1, the estimate θ k is kept fixed and the sufficient statistics S k 1 , S k 2;j and S k 3;j are computed sequentially using the current estimate θ k , the observations Y k and τ kþ 1 as the number of observations.The superscript k indicates which observations are used in the definition of the statistics.The next parameter estimate θ k þ 1 is computed when the last observation Y T k þ 1 is received using Algorithm 1.
Unlike in the traditional EM algorithm where the sufficient statistics are computed using forward-backward techniques, Cappé [5], Del Moral et al. [12] and Le Corff and Fort [22] proposed to compute the sufficient statistics recursively (i.e. as the observations are received and without any storage).In general state-space hidden Markov models, these online computations are not available in closed form (except in simple models such as linear Gaussian models) and have to be approximated, e.g. using sequential Monte Carlo methods [4,12].For finite state-space hidden Markov models, the computations can be performed in theory but are computationally too expensive if the number of states is too large (which is the case in our localization framework, see Section 5).Therefore, sequential Monte Carlo methods are used in this paper to localize the mobile.
In this case, the filtering distribution ϕ t k on the block k, i.e the distribution of where For all k Z 0 and t A f1; …; is a set of possible mobile positions at time T k þt.At each time step, the new population of particles is built from the previous population using the bootstrap filter, see Gordon et al. [20].The Bootstrap filter combines sequential importance sampling and resampling steps to produce this set of random particles with associated importance weights.Implementations of such procedures are detailed in Cappé [3], Del Moral [11], Cappé et al. [7], and Doucet and Johansen [14].
Online map estimation: We describe here the online approximation on the block Y k of the statistic S k 1 which is used to compute the map estimate θ k þ 1 when T k þ 1 is received.The computations for S k 2;j and S k 3;j follow the same lines.The rationale to establish this online approximation can be found in Del Moral et al. [12].The first particles ξ 0 i , i A f1; …; Ng, are sampled uniformly in C and the first weights are set as ω i 0 ¼ N À 1 , i A f1; …; Ng: Compute The approximation of S k 1 on the block Y k is then given by (iv) Once these computations are done for each statistic, the estimate θ k þ 1 is computed by Algorithm 1 applied with Ŝk 1 , Ŝk 2;j , Ŝk 3;j and τ k þ 1 .
The BOEM proposed in Le Corff and Fort [22] also introduced an averaged estimate f θk g k Z 0 based on a weighted mean of all the sufficient statistics computed in the past.It is proved in Le Corff and Fort [22] that this averaged estimator has an optimal rate of convergence.This can be easily done recursively after the computation of each statistic in step (iii).Step (iii) is then followed by ðiii 0 Þ Compute (with And step (iv) is then followed by ðiv 0 Þ Once these computations are done for each statistic, the estimate θk þ 1 is computed by Algorithm 1 applied with Sk 1 , Sk 2;j , Sk 3;j and T k þ 1 .Localization procedure: At each time step, we compute two estimators of the device position: (i) Nonaveraged algorithm: At each time step, the estimation of the device position is set as the particle with the greatest importance weight: ; where i max ¼ def argmax i ω i t : (ii) Averaged algorithm: As the average sequence f θk g k Z 0 has a better rate of convergence, another particle system fð ξi t ; ωi t Þg is run using only step (ii) (a) above by replacing θ k by θk .Then, the estimation of the device position is set as Other natural estimators of the device position such as the posterior means could be used: However, we did not observe significant differences between the estimated localization provided by the posterior means and by the proposed estimators with both simulated and true data.Moreover, we observed that the bootstrap filter sometimes produces several clouds of particles remote from each other (each of them is concentrated around local maxima of the filtering distribution).In this case, the proposed estimators offer a better accuracy than the weighted means.Finally, some indoor maps may not be convex such as in Section 5.2.In this case, the weighted means may not belong to the map while the proposed estimators always do.
Stabilization procedure: We add a stabilization step which consists in regularly replacing the original map estimate by the averaged one ðθ k ' θk Þ.This step is needed to improve the performance of the algorithm as detailed in Section 5.

Simulated data
In this section, all experiments are performed on the finite grid C ¼ f0; …; 30g Â f0; …; 30g.Note that in the following the particles are sampled on the square ½0; 30 Â ½0; 30 before being associated with the corresponding cells in the finite grid C. Each AP is modeled using the same coefficients c ⋆ and d ⋆ : Σ is a Gaussian covariance function defined by Σðx; Generalities about hyper-parameter determination in spatial data modeling can be found in Cressie [10].Ferris et al. [17] also describe the determination of the hyper-parameters v 1 and v 2 (when c ⋆ j and d ⋆ j are set to zero).In our case, these coefficients were calibrated after a measurement campaign in the office presented in Section 5.2.The corresponding grid stepsize is 1 m.Details about the calibration method for this particular model can be found in Dumont [15,Section 5.1].This measurement campaign might seem contradictory with our aim to get ride of such a campaign.However, we use this calibration to find relevant values for the true parameters c ⋆ j , d ⋆ j .We also use this measurement campaign to estimate the hyper-parameters v 1 and v 2 that characterize the prior distribution of δ ⋆ .Their values influence the smoothness of the Gaussian field δ ⋆ .An online estimation of v 1 and v 2 could be considered but is not addressed in this paper.We assess that a partial measurement campaign on a part of the environment only could be sufficient to calibrate them.We can also consider the same values of v 1 and v 2 for different environments so that v 1 ¼ 10 and v 2 ¼ 18 can be used directly and no measurement campaign is needed.The variance of the observation noise is s ⋆;2 ¼ 25 (this value was calibrated by computing the variance of a set of RSS observations at a fixed position).Thevariance of the transition kernel defined in ( 2) is chosen such that a¼6.
All runs are started with the same initial estimates: The number of particles is N¼ 25 and the initial position of each particle is chosen randomly and uniformly in C. The block sizes are given by Mapping error: For each map F j ¼ μ j þ δ j , the estimation error is set as the normalized L 1 error, such that the distance of a given map F j to the true map and the error displayed is the mean over all maps: Localization error: For a given block fT k þ 1; …; T k þ 1 g, the localization error is set as the 0.8-quantile of the sample: To assess the performance of our method we also display the estimated position given with a particle system run with the true maps F ⋆ j (i.e. using θ ⋆ instead of θ k in step ii) (a) of the algorithm).The localization error obtained using the true propagation maps F ⋆ and the true variance s ⋆;2 will be referred to as the reference estimator localization error.Fig. 4 displays the localization error for a different number of AP as a function of the number of blocks when the stabilization step is omitted.As expected, both the reference estimator localization error and the localization error of the averaged algorithm are improved as the number of AP increases.However, even with a great number of AP, the localization error of the averaged estimate does not converge to the reference estimator localization error.This is confirmed by Fig. 5 which displays the map estimation error and the localization error for the greatest number of AP (B ¼17).As shown in Fig. 5a the estimated position does not converge as the number of blocks (i.e. as the number of estimations) increases.After 50 blocks (about 40,000 observations) the position, which is badly estimated, does not provide good map estimates which increases the error on the averaged map estimate.Fig. 5b shows that both the map estimate and its averaged version do not converge.This convergence problem can be due to the curse of dimensionality since the number of parameters to estimate is high.Moreover, the higher the parameter space dimension is, the more likely the EM based algorithms are prone to converge towards local minima (see [7]).To overcome this difficulty, we propose to use the good behavior of the averaged estimate θk which offers a more accurate positioning than the nonaveraged version θ k (c.f.Fig. 5).Then θ k is regularly replaced by the averaged version θk .In Fig. 6, this stabilization process is performed each time N b ¼ 5 blocks have been used.As shown by Fig. 6a and b, this greatly improves the performance of the estimation of the maps and of the localization.Hence, the proposed algorithm uses this stabilization procedure and the averaged estimate to localize the mobile.
Figs. 7 and 8 illustrate the performance of the algorithm for the localization and for the estimation of the maps over 50 independent Monte Carlo runs.In Fig. 7, the reference localization error (i.e. when the maps are known) is also displayed.The convergence of the localization error to the reference error is almost reached after 100 blocks (about 100,000 observations).Similarly, the error for the estimation of the maps given by the averaged algorithm goes on decreasing after 100 blocks (the decrease is slower after 75 blocks).

True data
In this section, 10 AP are set up in an office environment.Fig. 9 represents a map of this environment as well as the position of the AP.The map is discretized using a grid C & f0; …; 32g Â f0; …; 28g.Note a major difference between the model given in Section 3 and the real data situation.For any measure Y k sent by the device, only several AP are represented in Y k .Therefore, the maps F j ¼ μj þ δj are not estimated simultaneously as, for any time step k, two AP might appear a different number of times in Y 1:k .We thus slightly modify our algorithm by introducing specific blocks and measure counters relative  to each AP.We shortly denote by "j A Y k " the fact that AP j belongs to Y k .
The variance s ⋆;2 is assumed to be known and its value ðs ⋆;2 ¼ 25 dBm 2 Þ is calibrated using a measurement campaign at a fixed position.About T ¼20,000 measures of the RSSI have been made by walking in the office for around 2 h and 45 min with a WiFi connected device (the device measures the RSSI every 0.5 s).Our algorithm produces position estimates however, unlike in the simulated data case, we do not have a direct access to the real position and thus cannot observe the localization error.To overcome this difficulty, we proceed to four phases of test.
During each phase, we regularly identify the true position associated with the last measure.For i in f1; 2; 3; 4g, we denote by S test;i the set of all the time steps belonging to phase i, and by fX t ; Y t g t A S test;i the data collected during this phase of test.These data will be used to compare the estimated positions f X t g t A S test;i with the true positions fX t g t A S test;i .We will also use these data to estimate the mapping error by considering, for any j ¼ f1; …; 10g, where 1 j A Yt equals 1 if measure Y t does contain AP j and equals 0 otherwise.Finally, we set We run our algorithm twice on data using the hyperparameters v 1 ¼ 10, v 2 ¼ 18.For these two experiments we will start the algorithm using different initial propagation maps: c 0 and d 0 being common to all AP and δ 0 being set to zero.
In the experiment 1, we consider initial parameters c 0 ¼ À37 and d 0 ¼ À9 which allow us to start the algorithm with initial estimators not far from the real propagation maps (see Table 1).In the experiment 2, we choose c 0 ¼ À37 and d 0 ¼ 9.In this case, d 0 being positive, the initial estimators state that the further the device is from an AP, the stronger the signal will be expected.This implies that the experiment 2 starts the algorithm with a completely wrong idea about how WiFi signals propagate in the environment (see Table 1).
In the two experiments the maps F j ¼ μj þ δj were updated a different number of times depending on the AP.This number of updates varies from two times for AP 2 to six times for AP 9.The evolution of the localization precision for the experiment 1 (resp.experiment 2) can be observed in Fig. 10 (resp.Fig. 11).Figs. 10 and 11 show that the localization improves with time for both experiments.We can observe in Fig. 10 that after a first period of improvement, the precision seems to reach a bound.On the contrary, with the experiment 2, the precision starts  really badly as we expected considering the initial point of our algorithm.The precision seems to stay constant until enough measures have been gathered by the device and until enough updates of the maps have been done.While the precision improves by around 1.5 m between the first test phase and the last one for the experiment 1, for the experiment 2, the precision considerably improves with time with a difference of 6.6 m between the beginning and the end of the experiment (see Table 2).Fig. 12 and Table 2 show that the experiment 2 final precision accuracy reaches (and even slightly overtake) the precision obtained with the experiment 1.
These observations confirm the robustness of our approach.If changes occur in the environment (modifying the way WiFi signals propagate), the difference between the resulting propagation maps and the current estimator will never be as substantial as it was between the true propagation maps and the initial estimate considered in the experiment 2. We reckon that our algorithm can adapt to such changes in order to the accuracy of the localization if the sufficient statistics S, k and T k are regularly reset.Finally, more elaborate localization algorithms could be considered to improve the accuracy.However, the aim of this paper lies in the efficiency of our propagation maps estimation procedure for localization purpose rather than the localization algorithm itself.

Conclusion
In this paper, we propose an online EM based algorithm to estimate the propagation maps needed in any WiFi based localization system.The main difference with the existing localization solutions is that these propagation maps are estimated using the data sent by the mobile device originally used for localization purposes.The existing WiFi based localization systems establish these propagation maps either in a deterministic way or by running a previous hand-made survey.In case of environmental modifications, the propagation maps are changed.Our technique can easily adapt to these changes by regularly reinitializing the estimation procedure while hand-made survey based systems cannot take into account these modifications without renewing the survey.Other elements could be analyzed such as the size of the environment or the materials constituting the obstacles in the environment that might particularly influence the "right choice" of the hyper-parameters v 1 and v 2 .An online estimation procedure of these hyper-parameters could be considered.It is possible to use different localization techniques using more sophisticated human motion model for instance in order to improve positioning accuracy.However, our algorithm only needs few information on the environment, namely the size of the indoor map and the AP positions.More sophisticated models might need  more information and thus make the installation process more complex.Table 2 Mean and 0.8-quantile (0.8-q) of the localization errors (in meter).

Fig. 1 .
Fig. 1.Histograms of the RSS frequencies for two AP.

Fig. 4 .
Fig. 4. 0.8-quantile of the distance between the true localization and the estimated position.The localization error is given with the averaged estimate (dotted line) and the reference estimate (bold line): (a) 5 AP; (b) 10 AP; and (c) 17 AP.

Fig. 5 .
Fig.5.Map estimation errors and localization errors with no stabilization step: (a) 0.8-quantile of the distance between the true localization and the estimated position.The localization error is given with the nonaveraged estimate (dotted line), the averaged estimate (dashed line) and the reference estimate (bold line) and (b) mean L 1 error on the map estimate with the nonaveraged estimate (dotted line) and the averaged estimate (dashed line).

Fig. 6 .
Fig.6.Map estimation errors and localization errors with the stabilized algorithm: (a) 0.8-quantile of the distance between the true localization and the estimated position with the stabilization process.The localization error is given with the nonaveraged estimate (dotted line), the averaged estimate (dashed line) and the reference estimate (bold line) and (b) mean L 1 error on the map estimate with the nonaveraged estimate (dotted line) and the averaged estimate (dashed line) with the stabilization process.

Fig. 7 .Fig. 8 .
Fig. 7. Boxplots of the localization error given by the stabilized algorithm with the averaged estimate (left) and the reference estimate (right) as a function of the number of blocks.

Fig. 9 .
Fig. 9. Map of the indoor environment with the position of each AP (dots) and their associated identification numbers.

Fig. 11 .
Fig. 11.Evolution of the localization precision for experiment 2: (a) tests errors and (b) tests errors CDF.

Fig. 12 .
Fig. 12. CDF of the last test phase errors for the two experiments.