Received signal strength–based indoor localization using a robust interacting multiple model–extended Kalman filter algorithm

Due to the vast increase in location-based services, currently there exists an actual need of robust and reliable indoor localization solutions. Received signal strength localization is widely used due to its simplicity and availability in most mobile devices. The received signal strength channel model is defined by the propagation losses and the shadow fading. In real-life applications, these parameters might vary over time because of changes in the environment. Thus, to obtain a reliable localization solution, they have to be sequentially estimated. In this article, the problem of tracking a mobile node by received signal strength measurements is addressed, simultaneously estimating the model parameters. Particularly, a two-slope path loss model is assumed for the received signal strength observations, which provides a more realistic representation of the propagation channel. The proposed methodology considers a parallel interacting multiple model–based architecture for distance estimation, which is coupled with the on-line estimation of the model parameters and the final position determination via Kalman filtering. Numerical simulation results in realistic scenarios are provided to support the theoretical discussion and to show the enhanced performance of the new robust indoor localization approach. Additionally, experimental results using real data are reported to validate the technique.


Introduction
The need for localization is not just confined to people or vehicles in outdoor environments, where Global Navigation Satellite System (GNSS) plays an important role for this purpose and is recognized to be the legacy solution. However, accurately estimating location indoors relying only on GNSS signals remains a difficult problem, mainly because of signal blockage or severe attenuations. Because of the advent of locationbased services (LBS) and the multiple technologies available for indoor solutions, there exists an actual need of robust and reliable indoor localization methodologies. 1,2 Due to the ubiquitous availability of powerful mobile computing devices, the boom of personalized context-and localization-aware applications has become an active field of research. A way of localization in indoor environments is using available radio signals such as wireless local area network (WLAN) (IEEE 802.11x), Zigbee, and ultra-wideband. The advantage of working with signals of the IEEE 802. 11 as the primary source of information to approach the localization problem is the inexpensive hardware and the already dense deployment of WLAN access points (APs) in urban areas. There are several channel models in the literature to characterize the indoor propagation environment. [3][4][5][6] In this article, the IEEE 802.11x model is considered because it does not require an accurate floor plan of the indoor scenario and can be implemented without using a third-party software. 7 Empirical and analytical models show that received signal strength (RSS) decreases logarithmically with distance for both indoor and outdoor channels. [8][9][10] This is the basis of the popular path loss model for RSS measurements. This model related the distance between a mobile node and the corresponding AP with the received RSS, being parameterized by the path loss exponent (quantifying the RSS attenuation with respect to distance) and the variance of a random term (modeling shadowing effects). In this article, we consider a two-slope path loss model 8,11,12 which extends the classical path loss model to account for the increase in path loss exponent and variance for large distances. The two-slope model can be described by two models which depend on a breakpoint distance value. 11 That is, for distances below the breakpoint, the RSS measurement obeys a first model with a given path loss exponent and variance, and above such breakpoint distance, a second model with its respective parameters. Compared to the path loss model, the two-slope model is able to better capture the complexity of signal propagation at large distances (i.e. above the breakpoint distance) to the AP, thus being more realistic.
In the literature, a plethora of different approaches and methodologies for geolocation and tracking a mobile node is available. 1,13,14 Typically, these works assume the one-slope path loss model where the parameters are assumed known, the latter being a strict assumption in real-world applications where the indoor channel model parameters are unknown to a certain extent. Interesting contributions proposed to use a twomodel approach in line-of-sight/non-line-of-sight (LOS/ NLOS) scenarios for the case of time-of-arrival (TOA)based localization, [15][16][17][18][19] and recently using RSS measurements for geolocation (static user) in outdoor urban scenarios. 20,21 Here, we are interested in algorithms that use RSS observations for locating and tracking the mobile node since most of mobile devices are equipped with wireless capability. 12 Indoor RSS-based localization has become a popular solution, but standard techniques still consider a time-invariant signal model with a priori known constant parameters. This standard RSS-based localization problem with known AP positions, a single-slope path loss model and known model parameters, has already been addressed in the literature using data fusion solutions. 13,22,23 While some contributions considered the RSS-based localization problem using a single path loss model with unknown parameters, 5,21,[24][25][26][27] the general solution considering the two-slope channel model is an important missing point.
While the single path loss model is adequate in free space propagation, a multi-slope piece-wise linear propagation model appears more suitable in indoor environments and in the presence of strong reflections. 8 This contribution considers the two-slope channel model and proposes a robust indoor localization solution based on a parallel architecture using a set of interacting multiple models (IMMs), 28-31 each one involving two extended Kalman filters (EKF) and dealing with the estimation of the distance to a given AP. Within each IMM, the two-slope path loss model parameters are sequentially estimated to provide a robust solution. Finally, the set of distance estimates is fused into a standard EKF-based solution to mobile target tracking.
The benchmarks to evaluate the performance of our IMM-EKF algorithm are the Crame´r Rao Lower Bound (CRLB) and the Posterior Crame´r Rao Lower Bound (PCRLB) derived in this work to provide a guidance in the improvement of the experimental design. The CRLB is used to assess the estimation of model parameters and the PCRLB the tracking solution. This, combined with a path loss exponent estimation, is a remarkable extension of the preliminary results presented in the work by Castro-Arvizu and colleagues. 32,33 The contribution of this article is a completely online two-slope channel calibration and, simultaneously, a mobile target tracking algorithm. The performance of the method is assessed through realistic computer simulations and validated real RSS measurements obtained from an experimental test.
The remainder of the article is organized as follows. The mathematical formulation of the system is given in section ''System model,'' and the problem formulation and main contribution are given in section ''Robust indoor localization: the problem.'' The proposed technical solution is detailed in section ''Adaptive IMMbased robust indoor localization solution.'' Illustrative simulations and results with real data are discussed in section ''Results,'' and section ''Conclusion'' concludes the article with final remarks.

System model
As already mentioned, the ultimate goal is the sequential localization (i.e. tracking) of a mobile device using RSS measurements from a set of N APs at known positions. A two-slope RSS model is considered to overcome the practical limitations of standard path loss models, which implies that the sequential position determination cannot be directly solved using a traditional filtering solution, and thus, a clever divide and conquer strategy must be adopted.
The proposed position estimation is performed in two steps: (1) estimation of relative distances to the set of visible APs; (2) fusion of these distance measurements into a blended tracking solution. This two-step estimation approach motivates the following general formulation. In this section, the peculiarities of the twoslope RSS model are presented, together with the statespace formulation for the distance estimation problem using RSS measurements and the location determination using distance measures, respectively.

Two-slope RSS measurement model
The widely used model for RSS observations is the path loss model, which is a simple yet realistic model for such measurements. It is parameterized by the path loss exponent (related to the power decay with respect to distance) and the shadowing (i.e. the random propagation effects). However, it has been observed in experimental campaigns that these parameters fluctuate and are indeed distance dependent. As a conclusion, the parameters employed in the traditional path loss model are highly site-specific. 34,35 Therefore, in many situations, more sophisticated models should be required.
In this work, an extension of the classical path loss model accounting for two regions of propagation, referred to as the two-slope model, 11 is considered to overcome the practical limitations of the standard case. Path loss is the reduction in signal strength over distance. The path loss depends specifically on the distance between the transmitter (Anchor Point) and the receiver (Mobile Target). 36 Indeed, it has been observed that for far distances (5 d 30 m), there exists a steeper overall drop in the RSS at the receiver. This effect is due to path reflections from the environment (specially reflections from walls and ceilings).
Under this model, the RSS for the rth AP (r = 1, 2, . . . , N total number of APs) to the mobile target is modeled as 11 where d is the relative distance between the AP and the moving node where the RSS was measured, and h (2) where L 0 is the RSS in a reference distance. 11,37 The first equation gives the path loss (in dB) for close distances (d d bp , known as the breakpoint distance), and the second equation gives the path loss beyond d bp . The a 1 and a 2 values are referred to as path loss exponents, defining the slopes before and after d bp , respectively. In the sequel, we refer to these two RSS models as M 1 and M 2 , respectively. The functions h (1) (d) and h (2) (d) were obtained by measurement campaigns using radio signal ray tracing methods, premeasured RSS contours centered at the receiver, or multiple measurements at several base stations. 12,22,38 Depending on the transmitter/receiver geometrical configuration, the RSS measurements might be distorted from the nominal. This variation (known as shadow fading or log-normal shadowing) can be modeled by an additive zero-mean Gaussian random variable. The notation x s 2 ;N (0, s 2 ) is used. As it happens for the path loss exponents, the variance values differ before and after the breakpoint distance. The standard deviations of the received power before and after breakpoint distance, s 1 and s 2 , are expressed in units of dB and are assumed relatively constant with distance. To sum up, the two-slope RSS measurement model is parameterized and fully determined by the following set c d = ½a 1 , a 2 , s 2 1 , s 2 2 , d bp T .
State-space formulation: RSS ! distance As previously stated, the proposed strategy to solve the localization problem uses a two-step approach. In the first step (i.e. distance estimation) and for the rth AP, the observations correspond to the RSS measurements and the unknown states to be sequentially inferred are where d r k is the distance between the mobile and the rth AP and _ d r k is the rate of change in this distance. The linear evolution of states is u r k = A u u r kÀ1 + v r k where v r k is the process noise accounting for possible modeling mismatches, such as a possible acceleration of the mobile. In other words, this noise term gathers different forces that could affect target's dynamics and which are not explicitly modeled. The process noise is normally distributed with zero mean and covariance matrix Q k = s 2 d B u B T u 24 where s 2 d models the uncertainty on the mobile dynamics. The state equation includes these matrices where Dt is the sampling period.
To complete the state-space representation, the observation vector is defined. In this case, the RSS measurements per AP are precisely the observations used to infer u r k , and thus, y r where we recall that the model for RSS r (d) depends on the breakpoint distance. Therefore, h(d k ) has to be selected according to equation (1), and variance of the measurement noise n k is s 2 1 or s 2 2 . To conclude, the state-space formulation for the pair fRSS, dg and one single AP is given by u r k and y r k . If the state estimation or filtering problem taking into account this state-space formulation is to be solved using an EKF solution, the following Jacobian matrices of the measurement functions are needed because h (1) and h (2) are nonlinear. Where the corresponding 2 3 1 Jacobian matrices are as follows The final goal is to sequentially obtain the mobile position. The location calculation is performed using the N distance measures related to the N visible APs. In this case, the state vector gathers the mobile position and velocity, x k = ½x k , y k , _ x k , _ y k T , and the observation vector is z k = ½d 1 k , . . . ,d N k T ; in this case,d r k refers to the noisy distance measurements between the target and the rth AP at time k. The state-space mathematical representation of the problem is given in the sequel.
Process equation. The state evolution is given by p is the variance related to the mobile acceleration, and Observation equation. The relation between distance and position is defined as z k = g k (x k ) + n k , where z k 2 R N is a vector gathering the estimated distances to the N APs, the observation error n k is modeled as an uncorrelated white Gaussian noise with covariance R p, k , and the nonlinear observation function, g k (x k ), is defined as the distance of the mobile to every anchor point where fx r p , y r p g is the position of the rth AP, d r k is the true distance between the mobile and the rth AP, and the corresponding Jacobian used to solve the nonlinear filtering problem into an EKF solution is given by Robust indoor localization: the problem In the previous section, the overall system model has been fully described, providing the mathematical formulation needed to come up with a new powerful and robust indoor localization solution. Notice that the first state-space model is parameterized and completely determined by , d bp T , and the second one, by c (2) = ½s 2 p , R p, k T , where the positions of the N APs are assumed to be known. To understand the problem at hand and to clearly place the new contribution of this work with respect to stateof-the-art RSS-based localization solutions, the main estimation problems within the indoor localization framework are summarized in the sequel: Generalized localization problem. In this contribution, the term ''generalized'' refers to the use of a generalized two-slope path loss model (still considering perfectly known model parameters) to overcome the practical limitations of the conventional approach. This can be seen as an improvement in the standard localization case. To solve this problem, a more sophisticated solution is needed to cope with the different distance-dependent RSS measurement models. To the best of our knowledge, no contributions in the literature face this generalized localization problem; thus, the solution provided in this work can be simplified to optimally cope with it when model parameters are known. Robust localization problem. In practice, the system parameters are not perfectly known, what leads to poor localization performances when using the previous solution in real-life applications. The robust or adaptive localization problem implies the sequential position determination from a set of RSS measurements and the simultaneous model calibration or system model parameter estimation. Two cases can be considered: -Case I: standard model. Few works have considered the localization problem using RSS measures together with the adaptive estimation of the propagation model parameters, 5,21,[24][25][26][27] and these contributions only take into account simple propagation models. -Case II: generalized two-slope model. A general solution to the robust RSS-based indoor localization problem considering a generalized realistic propagation model is an important missing point in the literature because the position determination assuming known system parameters and/or simplified system models is not of practical interest. This is the main problem of interest in this contribution.
This article proposes a solution to the generalized robust indoor localization problem, considering a twoslope RSS propagation model and the estimation of the path loss exponents and shadowing of the model. Mathematically, this is expressed as the sequential esti- to N, and the simultaneous estimation of ½a 1 , a 2 , s 2 1 , s 2 2 T as a model calibration strategy. Why the other parameters, are not estimated and how they are determined will be discussed in the following sections.
Taking into account the robust filtering problem at hand, the general system model formulation, and the two-step estimation approach, the following points are sequentially treated in the sequel: (1) distance estimation from a set of RSS measures (section ''Parallel IMM-based solution for distance estimation''), (2) position estimation from the intermediate distance estimates (section ''Position determination''), and (3) model calibration and parameter estimation strategies (section ''Maximum likelihood model calibration''). To close the loop, the overall discussion on the complete robust localization solution is given in section ''Overall robust IMM-based architecture with ML-based system parameter estimation.''

Adaptive IMM-based robust indoor localization solution
The main concern of this section is to provide a clear answer and the mathematical derivation and solution to the problem stated in the previous section. The first goal is to justify the reasoning behind the use of an IMM-based approach to solve the distance estimation problem from a set of N RSS measurements. Then, to detail the proposed solution to obtain the final position from these intermediate distance estimates. Together with the main architectural core, the model calibration strategies and the corresponding maximum likelihood estimators (MLEs) for the two-slope model parameters are derived. To sum up, at the end of the section, the overall architecture is given.

Parallel IMM-based solution for distance estimation
The first approach that comes to mind to solve the distance determination problem is the use of a traditional filtering solution, such as the EKF, where the observation accounts for the full set of RSS measurements y k = ½y 1 k , . . . , y N k T and the global state evolution takes into account the N individual states, u k = ½u 1 k ; . . . ; u N k . But it is straightforward to see from equation (1) that this is not a valid approach because the measurement model is distance dependent. The two-slope model can be used to effectively model the distance between the AP and the mobile node but both implicit models must be treated separately. The natural solution to overcome this model-switching problem is to use an IMM-based approach.
The key idea behind the IMM is to use a bank of L KFs, each one designed to cope with a specific model (or model set), and to obtain the state estimation as a clever combination of the individual estimates. If the full set of N independent observations y k is considered, the question that arises is how many KFs should be considered into the IMM as each independent observation may obey M 1 or M 2 , the answer is 2 N filters (i.e. all the possible combinations of M 1 and M 2 for the N observations). It is clear that this is not a practical solution for an arbitrary number of APs; therefore, again a divide-and-conquer strategy treating independent measurements separately is the best solution. In this contribution, a parallel IMM-based approach is adopted, considering N IMMs each one involving two KFs according to the two path loss models. The block diagram of a single IMM is sketched in Figure 1.
At each discrete-time instant k, the IMM algorithm follows a clear three-step architecture: (1) Reinitialization, which represents the interaction between filters; (2) KF, being the individual filtering solutions; and (3) Model probability, which computes the model likelihood to decide whether the input measurement comes from M 1 or M 2 . The final estimates are obtained as a weighted combination of the individual KF outputs using the corresponding model likelihoods. Mathematically, one cycle of the standard IMM associated with the rth RSS measurement is sketched in Algorithm 1, where p ji (for i, j = 1, 2) is a two-state Markov model transition probability matrix to describe the switching system. The Markov chain is shown in Figure 2.
Notice that in its standard form, both the two-slope model parameters and the process noise variance, gathered in vector c (1) , must be specified in the IMM. These parameters must be set to the true ones for an optimal solution. Moreover, the initialization of both EKFs and each AP, fû (i), r 0j0 , P (i), r 0j0 for i = 1, 2g, must be set according to the problem at hand. The following section is an example of a possible parameterization (i.e. these values are the ones used later in the simulations).
The adaptive estimation of the two-slope model parameters to obtain a robust filtering solution is provided in section ''Maximum likelihood model calibration,'' together with the discussion on how to set s 2 d and d bp . The error covariance matrix has an initial value assigned as P (i), r 0j0 = 4Q k for each AP. The initial value state vector for the filter isû (i), r 0j0 = u r 0 + v where v;N (0, 0:8I).

Position determination
At each time step k, the output of the rth IMM filter is an estimate of the distance (and distance rate) to the rth AP,û r kjk . Therefore, the output of the whole block of parallel IMMs is an estimated distance vector, z k = ½d 1 k , . . . ,d N k T , whered r k is the first element inû r kjk . Using the state-space formulation given by x k and z k , and taking the estimated distance vector as input measurements, the position determination can be solved straightforwardly with an EKF.
The EKF-based position estimation is sketched in Algorithm 2. Notice that the measurement Jacobian matrix in equation (7) is evaluated at the predicted state to obtain a linear formulation, thus, G k = G k (x k )j x k =x kjkÀ1 . Regarding the implementation of such solution, both the initialization, fx 0j0 , P x, 0j0 g, and the noise statistic parameters, fs 2 p , R p, k g, must be specified: The measurement noise covariance can be set according to the output of the individual IMM filters. The measurement noise accounts for possible errors or observation noise, but the set of measurements used in the position determination is the set of estimated distances. The error on the estimation of the distance within the  IMM is given by the estimation error covariance matrix, which for each AP is P r kjk . Therefore, at each time step, the measurement error covariance is set according to R p, k = diag (P 1 kjk (1, 1), :::, P N kjk (1, 1)).
The error covariance matrix has an initial value assigned as P x, 0j0 = 4s 2 p B x B T x for each AP. The initial value of s 2 p is indicated in section ''Maximum likelihood model calibration.'' The initial value state vector for the filter iŝ x 0j0 = x 0 + v, with v;N (0, 0:8I).

Maximum likelihood model calibration
In the previous paragraphs, the proposed two-step state estimation solution has been derived, first providing an intermediate distance estimate using a bank of N IMMs to cope with the N APs and the two-slope path model and then an EKF-based solution to obtain the final position from these estimates. Both state-space models are determined by a set of parameters, c (1) and c (2) , respectively, which in general may be unknown to a certain extent. It has already been justified that the model parameters must be somehow adjusted or estimated to obtain a robust and flexible solution. The second state- Algorithm 1 Step k of the IMM for the r th AP.
1: For i = f1, 2g and j = f1, 2g 2: Reinitialization: Calculation of the predicted mode probability, mixing weights, mixing estimates, and mixing covariances, respectively 3: Model-conditioned EKF: Prediction, innovations' covariance matrix, Kalman gain, state estimate, and the corresponding error covariance matrix are given byû 4: Model probability update: The model likelihood function and model probability are, respectively Algorithm 2 EKF formulation for position determination.

7: Estimate the updated statê
x kjk =x kjkÀ1 +K k z k Àẑ kjkÀ1 À Á : 8: Estimate the corresponding error covariance: P x, kjk = P x, kjkÀ1 ÀK k G k P x, kjkÀ1 : 9: Set k ( k + 1 and go to step 2. space model parameters' specification, c (2) , has already been detailed in section ''Position determination,'' and thus hereafter the focus is on the sequential estimation of c (1) . This section presents the estimation strategies for the two-slope model calibration together with the theoretical derivation of the ML (maximum likelihood) parameter estimators for ½a 1 , a 2 , s 2 1 , s 2 2 T . Recall that the first state-space model has two extra parameters, namely, the process noise variance s 2 d and the breakpoint distance d bp , which are not estimated but rather set as detailed in the sequel: The process noise variance for the distance estimation problem, s 2 d , is determined assuming that the mobile target has an average velocity of 1 m/ s in the simulations presented in this work, so a small initial value for s 2 d of 0.7 m/s 2 is chosen. The same reasoning applies for s 2 p . The breakpoint distance indicates the change in models in the path loss model. An off-line Bayesian approach for the d bp estimation was considered in Castro-Arvizu et al., 37 but its applicability to the on-line configuration of interest here is still under investigation. In this work, the IEEE standard 34 channel parameters are considered to set the value of d bp = 5 m. For the sake of completeness, in the computer simulations, the impact of under-and over-estimation of this value is presented.
In the two-slope model (1), the RSS measurements may come from the first equation modeling the propagation for close distances (called M 1 ), or alternatively, they may obey the second equation modeling the propagation for distances beyond the breakpoint distance (called M 2 ). Therefore, at the input of the system, there exists a model uncertainty which must be resolved.
In the proposed methodology, each IMM inherently treats this model uncertainty by computing the model likelihood from the innovations of each KF. For each AP r and model i, the model probability is given by h (i), r k . These probabilities are used in the filter to weigh the outputs of the individual KFs but can also be reused for the model calibration. At each time step and using these model probabilities, two subsets of RSS measurements are constructed: if h (1), r k .h (2), r k , the RSS measurement y r k is associated with y r 1, k (i.e. which represents the RSS measurement subset obeying M 1 ), Y r 1, k = fy r k : h (1), r k .h (2), r k g; otherwise, it is associated with Y r 2, k (i.e. which concatenates the RSS measurements obeying M 2 ). Therefore, we take a hard decision to associate RSS measurements with each of the models and construct the MLE accordingly. Note that the cardinality of these sets is upper bounded by the present time instant, U i ¼ D jY r i, k j\k, i = f1, 2g, and that their sum is precisely jY r 1, k j + jY r 2, k j = k. For the sake of clarity in the forthcoming derivations, let us define the elements in the sets as Y r 1, k = fy r 1, 1 , . . . , y r 1, U 1 g and Y r 2, k = fy r 2, 1 , . . . , y r 2, U 2 g: MLEs for M 1 . In the following, the MLEs for the close distance (d d bp ) model parameters are provided. Notice that the parameters to be estimated in this case are fa 1 , s 2 1 g at each AP. As the parameter a 1 appears in both models,â 1, 1 is used for the close distances and a 1, 2 for distances beyond d bp , both being estimators of a 1 .
The expressions for the estimators are herein provided, the reader is referred to Appendix 1 for the complete derivation. For the rth AP, the 'th sample of the first model subset Y r 1, k at time k is Gaussian distributed, y r 1, ' ;N ( y r 1, ' , s 2, r 1 ), with y r 1, ' = L 0 + 10a r 1, 1 log 10 d r ' . Using this subset and assuming a known distance to the rth AP, d r ' , at instant k, the MLEs are given bŷ a r 1, 1 = MLEs for M 2 . Following the same procedure but using the second model subset of RSS measurements, the MLEs of ½a 1, 2 , a 2 , s 2 2 are given in the sequel for the rth AP. Refer to Appendix 2 for the complete derivation.
For the rth AP, the 'th sample of the second model subset Y r 2, k at time k is Gaussian distributed as well, y r 2, ' ;N ( y r 2, ' , s 2, r 2 ), with y r 2, ' = L 0 + 10a r 1, 2 log 10 d bp + 10a r 2 log 10 (d r ' =d bp ). Using this subset and assuming a known distance to the rth AP, d r ' , at time instant k, the MLE of s 2 2 is given byŝ and the MLEsâ r 2 andâ r 1, 2 are the first and second elements of vectorâ in equation (33).
For the MLEs of both models, there are two issues to account for. First, notice that in practice, the true distance to the rth AP is unknown, and thus, an estimated r ' must be used instead, which is available at the output of the corresponding IMM.
Second, two ML estimates exist for a r 1 , but the algorithm needs to take a decision and provide a unique solutionâ r 1 . In this work, the following simple rule is considered. The idea is to use the estimator that has Overall robust IMM-based architecture with ML-based system parameter estimation This subsection summarizes the overall system. A diagram of the interconnection of elements can be consulted in Figure 3.
The main goal of the proposed method is to use the proposed IMM-EKF algorithm to sequentially estimate both the two-slope path loss model parameters and the mobile target position. Notice that in practice, these parameters are typically found after a scene analysis and the posterior linear regression on a semilogarithmic scale. 3,10,34,35,37,39 To avoid this off-line site-dependent procedure, the proposed solution performs the on-line parameter estimation within the IMM architecture.
In general, the IMM is capable of dealing with model transition, which are modeled as a two-state Markov jump process, being used to filter the distance measurements. The state u r k is simultaneously estimated by the two EKFs in the IMM according to the corresponding model, which are fused using the model probabilities to obtain the final estimate. As already mentioned in section ''Maximum likelihood model calibration,'' the model probabilities h (i), r k resolve whether the RSS measurement acquired at instant k obeys M 1 or M 2 . These probabilities are used to construct the two measurement subsets to obtain the ML path loss parameter estimation, which in turn are feed back to the corresponding EKF. In the proposed architecture, an IMM-EKF is run for every AP. The set of distance estimates is used to form the blended positioning solution. The EKF used to obtain the final mobile position was described in section ''Position determination.''

Evaluation of the IMM-based robust indoor localization algorithm with real data
This section evaluates the proposed algorithms (calibration and distance filtering) with real RSS measurements.

Experimental setup
The measurements were obtained in a real office environment as shown in Figure 14 in a NLOS scenario. The architectural plan is the second floor of a typical multi-story office building with drywall and wood wall panelings reinforced with aluminum bars. These RSS measurements were used for the model calibration and distance estimation also. For the tracking task, two algorithms were employed: an algorithm using an EKF modeled with the classical one-slope model and another with our proposed IMM-EKF algorithm.
The mobile path is a test-bed used to collect the RSS measurements that includes a RaspberryPi board. The test-bed and the whole hardware description are presented in section ''Hardware description.''

Hardware description
The ranging/positioning payload is a development board with multiple connections where ranging and positioning algorithms can be easily implemented. The RaspberryPi board is the model B with a Universal Serial Bus (USB) WiFi card (IEEE 802.11n, 802.11g, 802.11b). The RaspberryPi board has a Central Processing Unit (CPU) ARM11 @700 MHz featuring with a floating point Arithmetic Logic Unit (ALU), Ethernet, ALU 2.0, I2C bus, a serial port and generalpurpose input/outputs (GPIOs), a Linux Operative System (debian-based distribution), and a dedicated high-definition camera connector. The positioning payload is a cheap and easy-to-use system that allows WiFi RSS reading with a CPU power (similar to an entry level smartphone).
The overall system consists of the ranging/positioning payload and the data base. The RaspberryPi microcontroller sends the RSS measurements to the data base of the server. In Figure 4, a schematic of the overall system is presented where the ranging/positioning payload reads RSS WiFi measurements employing a TL-WN722N WiFi card (from TP-LINK manufacturer).
An integrated navigation information system must continuously know the current position with a good precision, and thus, a model is needed to measure the real position. The chosen model is a four-wheel Robot that is capable of performing a programmed trajectory through waypoints. The main board features an Arduino Platform. Arduino is an open-source platform and consists of a physical programmable circuit board (often referred to as a microcontroller) and an IDE (Integrated Development Environment) based in C+ + language programming and used for software loading in the board. The initial value forû 0j0 is the first real distance measurement given by the Robot.

Results
In section ''Simulation results,'' the results of the IMM-EKF algorithm were obtained with synthetic signal and in section ''Validation with real data'' with real data.

Simulation results
The method proposed in this work was validated by computer simulations in a scenario depicted in Figure 5 that could be considered as a realistic scenario and where the number of APs (N = 6) were deployed in a 30 3 30 m 2 area at known locations. A mobile node was moving in the area, with initial coordinates being at origin, (0, 0), and a steady velocity of 1 m/s corresponding to a pedestrian movement. The duration of the trajectory consisted in 18 s with a sample period Dt of 100 ms, emulating realistic WLAN measurements. The trajectory was kept the same (see Figure 5) throughout the simulations in order to see particular cases such as node approaching an AP (e.g. AP 2), leaving the proximity of an AP (e.g. AP 5), and combinations (e.g. AP 6). These three APs (e.g. 2, 5, and 6) are those where model switching could take place. The rest of them have a relative distance to the mobile node larger than d bp = 5 m along the track, and thus, they are likely to provide observations for M 2 . The true model parameters were a 1 = 2, s 1 = 3, a 2 = 3:5, and s 2 = 5.
First, a batch of simulations for a single realization of the method is shown which allows us to comment details, as well as to provide some insights and intuition regarding the operation of the proposed method. Second, Monte Carlo simulations were performed to evaluate the root mean square error (RMSE) performance of our method, compared to other comparative solutions. Namely, we compared our method (termed in the legends as IMM-MLE) with the same IMM architecture able to track the node under the two switching models but with known model parameters, in which case the MLEs are not required since the true values were set. Also, we compared the solution to a standard EKF-based algorithm considering that all observations obey M 1 and thus not accounting for the possible changes in the path loss exponent and shadowing variance for relative distances larger than d bp . The latter method assumed perfect knowledge of M 1 parameters.
For a single realization, the distances estimated compared with the real ones with respect to every AP along the simulation duration are shown in Figure 6. For better understanding of these results, the d bp value is also plotted. Notice that in some time instants, the distance between mobile target and some APs changes from values above d bp to lower values and vice versa. This indicates that a switch model has occurred. The IMM structure is able to track the distances even under model switching.
In the situation where the mobile node is close to the breakpoint distance (i.e. the border for the two-slope model), the model probabilities h (1), r k and h (2), r k are likely to exhibit nervous behaviors. To illustrate this fact, the estimated distance referred to AP 5 along the simulation is shown in Figure 7. The bottom plot therein shows the performance of the decision process in M 1 -M 2 switching. The top plot presents the estimated distance. Figures 8 and 9 show the MLE of the model parameters for AP 5, featuring their convergence to the true values. For the sake of completeness, the bottom plot corresponds to the computed model probabilities. The decision process between updating the estimatorâ 1 (either according toâ 1, 1 orâ 1, 2 ) is visible in the elapsed time when the model target remains after the model switching.
The RMSE performance was evaluated and compared to the benchmark methods detailed earlier.    Figure 10 shows the average RMSE of distance estimation over all six APs. Figure 11 shows the corresponding RMSE and PCRLB of position estimation, the latter computed recursively as in Tichavsky et al. 40 From these figures, it is highlighted that our robust indoor localization method has good accuracy when compared to a method that has full knowledge of the model parameters. Clearly, the standard method (optimally designed to operate on M 1 ) cannot handle model mismatch on RSS measures. This is mostly because there is a large probability of being within M 2 than in M 1 . Focusing on AP 5 as our illustrative link, the corresponding RMSE of the four parameter estimation is shown in Figure 12, where we observe the convergence of the ML-based method after some instants. Finally, we performed some sensitivity tests of the proposed method to deviations from the assumed model parameters. Namely, the designed robust     Refer to an explanation in section ''Maximum likelihood model calibration.'' However, this instance is taken as a sensitivity analysis of our algorithm to these values. Both cases were taken separately. The first analysis is when d bp = 3 m (underestimation) and d bp = 10 m (overestimated), recall that d bp = 5 m in the simulated data. Figure 13 shows the RMSE of position estimation under d bp mismatches, showing that the algorithm has more sensitivity when the value of d bp is overestimated. For the case of having s 2 d deviations of 10s 2 d and 100s 2 d , we provide the average RMSE over the simulated trajectory, for the sake of brevity. The average RMSE is 0:3571 m, 0:47551 m, and 0:5701 m, for true s 2 d , 10s 2 d , and 100s 2 d , respectively. These results confirm that the method is not very stressed when s 2 d is not perfectly known.

Validation with real data
The IMM-EKF algorithm was implemented to estimate the distance to the AP in an NLOS environment as shown in Figure 14. The error between the estimated distance and the real mobile target position per every distance interval is shown in Figure 15 for all deployed APs. To verify the accuracy of our algorithm, the distance estimation was also computed with an EKF modeled with the classical one-slope model only. From this figure, it is notable that for large distances, the IMM-EKF algorithm has a best performance than just considering the classical path loss model only.
The mobile path tracking is shown in Figure 14 and the distance estimation error to every AP can be seen in Figure 15. With real RSS measurements and an on-line channel parameterization, the implementation of the IMM-EKF algorithm has a good consistency in positioning terms that can be demonstrated in Figure 16.  For the sake of comparison, Figure 17 shows the error cumulative distribution function (CDF). Our EKF-IMM algorithm gives an improvement of 53.15% and a mean error of 0.1934 m in comparative with using only the EKF modeled with the classical oneslope model that gives a mean error of 0.4128 m.

Conclusion
Mobile location via RSS measurements has been formulated as a switching nonlinear state-space problem, accounting for realistic conditions where RSS measurements were seen to follow two propagation models depending on the relative distance to the reference nodes. This work proposed a robust IMM-EKF algorithm, including an on-line ML estimation procedure to sequentially adapt the model parameters. Using the model likelihood functions, the proposed method provides accurate distance estimates between the mobile and each anchor point, which are used for position tracking. Simulation results under realistic WLAN scenarios showed that the proposed IMM-EKF algorithm provides both a good mobile estimation and channel calibration, and much better performance compared to the state-of-the-art techniques. An analysis of the breakpoint distance sensitivity was made in this work but as future work, an on-line breakpoint distance estimation is proposed as well as an NLOS/LOS algorithm identifier.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been supported by the Spanish Ministry of Economy and Competitiveness through project TEC2015-69868-C2 Applying natural logarithm, the log-likelihood is log L(u 1 jY r 1, k ) = À U 1 2 log (2ps 2 1 ) À 1 The maximum likelihood estimator (MLE) of u r 1 is obtained from the following what leads to the following estimator Then, the path loss exponent, a 1, 1 , is computed from ∂ log L(u 1 jY r 1, k ) ∂a 1, 1 what leads to the following estimator