Abstract

This paper extends an algorithm that exploits multipath propagation for position estimation of mobile receivers named Channel-SLAM. Channel-SLAM treats multipath components (MPCs) as signals from virtual transmitters (VTs) and estimates the positions of the VTs simultaneously with the mobile receiver positions. For Channel-SLAM it is essential to obtain angle of arrival (AoA) measurements for each MPC in order to estimate the VT positions. In this paper, we propose a novel Channel-SLAM implementation based on particle filtering which fuses heading information of an inertial measurement unit (IMU) to omit AoA measurements and to improve the position accuracy. Interpreting all MPCs as signals originated from VTs, Channel-SLAM enables positioning also in non-line-of-sight situations. Furthermore, we propose a method to dynamically adapt the number of particles which significantly reduces the computational complexity. A posterior Cramér-Rao lower bound for Channel-SLAM is derived which incorporates the heading information of the inertial measurement unit (IMU). We evaluate the proposed algorithm based on measurements with a single fixed transmitter and a moving pedestrian carrying the receiver and the IMU. The evaluations show that accurate position estimation is possible without the knowledge of the physical transmitter position by exploiting MPCs and the heading information of an IMU.

1. Introduction

Today, most smartphones are equipped with global navigation satellite systems (GNSSs) receivers which allow using applications on the smartphones for navigation [1]. GNSSs provide sufficient position accuracies for mass market application in open sky conditions. However, indoors or in urban canyons the GNSS positioning accuracy could be drastically reduced. In these situations, the GNSS signals might be blocked, degraded by multipath effects, or received with low power. To enhance the positioning performance indoors, different methods and sensor systems can provide position information rather than relying on GNSSs [2–4]. Most of the indoor positioning systems use local infrastructure like positioning with Radio Frequency Identification (RFID) [5], mobile communication base-stations [6, 7], wireless local area network (WLAN) [8], or ultra-wideband (UWB) [9–11]. However, also these wireless radio technologies experience multipath and non-line-of-sight (NLoS) propagation. Multipath propagation is experienced when the transmitted signal arrives at the receiver via several propagation paths. These propagation paths with different delays are caused by reflections, diffractions, and scattering of the electromagnetic wave. Hence, the signal at the receiving antenna consists of a superposition of multiple replicas of the transmitted signal, where each version is called multipath component (MPC) traveling along an individual propagation path. The delay estimate of standard algorithms like the delay locked loop (DLL) is biased in multipath propagation environments [12]. Algorithms like [13–15] reduce the multipath error by modifying the DLL structure. Other algorithms estimate the channel impulse response (CIR) in order to mitigate the influence of multipath propagation on the delay estimate, for example, [16–20]. To retrieve the required delay from the CIR, the path with the smallest delay is treated as the line-of-sight (LoS) path. However, treating the smallest delay as the LoS path may result in weak positioning performance in NLoS situations. Furthermore, even advanced multipath mitigation algorithms reduce the multipath effects only to a certain degree due to limited signal bandwidth and measurement noise [18].

Nowadays, multipath exploitation instead of mitigation is attracting more and more interest. The authors of [21, 22] exploit multipath propagation for positioning of mobile terminals using multipath fingerprinting algorithms. Other algorithms, for example, [23, 24], interpret reflected signals as signals emitted from virtual transmitters (VTs), where the VT positions are precalculated based on the knowledge of the reflecting surface and physical transmitter positions. Furthermore, the authors of [25] estimate and track the phase information of MPCs using an extended Kalman filter (EKF) and estimate the user position using a time difference of arrival (TDOA) positioning approach. Other algorithms like [26] use a nonlinear least squares algorithm combining UWB measurements at several receiver positions to estimate the positions of the VTs and the receiver simultaneously within small scale scenarios.

This paper describes and extends the multipath assisted positioning algorithm referred to as Channel-SLAM; see [27–31]. Channel-SLAM considers a moving receiver and is suitable for GNSS denied areas like indoor areas. Similarly to other multipath assisted positioning approaches, Channel-SLAM interprets MPCs as LoS signals emitted from VTs. In addition to reflected signals, Channel-SLAM considers also paths occurring due to multiple number of reflections, diffractions, or scattering as well as combinations of these effects. As a consequence, the reception of several MPCs allows position estimation even if only one physical transmitter is present. Interpreting MPCs as directly propagated signals originated from VTs, Channel-SLAM enables positioning also in NLoS situations. Additionally, Channel-SLAM does not require any prior knowledge on locations of reflecting surfaces as Channel-SLAM estimates the receiver position, velocity, clock bias, and the VT positions simultaneously which can be interpreted as simultaneous localization and mapping (SLAM) with radio signals. In [27, 28, 31], we showed that positioning is possible in NLoS scenarios using MPCs without the knowledge of the room geometry by using Channel-SLAM. We investigated in [27] TDOA positioning and especially TDOA between MPCs such that time synchronization between physical transmitters is not essential. In [31], we derived Channel-SLAM based on a Rao-Blackwellized particle filter (RBPF) and compared the accuracy of Channel-SLAM to a derived posterior Cramér-Rao lower bound (PCRLB). However, the Channel-SLAM algorithms in [27, 28, 31] use linear antenna arrays and assume the knowledge of the physical transmitter position.

In this paper, we propose an implementation of Channel-SLAM that uses only a single receiving antenna and fuses similarly to [29, 30] additional information obtained from an inertial measurement unit (IMU). Today many smartphones feature Microelectromechanical System (MEMS) IMUs, which can provide short term relative orientation and position information. Theoretically, the measurements of the IMU can be directly used in an inertial navigation system. However, the position calculation involves double integrations; hence, even small measurement errors quickly cause a drift in the position solution [32]. To avoid that, we only fuse heading measurements from the IMU which solely requires an alignment of the coordinate systems. The heading information of the IMU allows improving the performance of Channel-SLAM by resolving ambiguities and angle of arrival (AoA) measurements are not mandatory anymore. Being a relative positioning system, Channel-SLAM requires an initial prior knowledge of the receiver position and moving direction to define the coordinate system. The positioning algorithm derived in this paper is based on a RBPF where we employ a new transition model for pedestrians. In [29, 30], we showed that positioning with only one physical transmitter is possible if MPCs and heading information from an IMU are used. Compared to [29, 30], the novel transition model enables a performance gain in the position accuracy. In addition to [27–31], we propose a method to dynamically adapt the number of particles which significantly reduces the computational complexity. Furthermore, a PCRLB for Channel-SLAM is derived which incorporates heading information obtained by using an IMU. The developed positioning algorithm is evaluated based on measurement data obtained in an outdoor scenario, where the position of the physical transmitter is unknown. Based on these measurements, we compare the accuracy of Channel-SLAM to that of the derived PCRLB.

The paper is structured as follows: Section 2 describes the signal model; afterwards, Section 3 describes the proposed algorithm which is split into four subsections: Section 3.1 addresses Channel-SLAM; Section 3.2 describes two different transition models using the heading information from an IMU; Section 3.3 summarizes the RBPF; Section 3.4 describes the implementation of the RBPF; afterwards, we derive in Section 4 the PCRLB for Channel-SLAM incorporating the heading changes of the IMU. Thereafter, Section 5 evaluates the algorithm based on measurement data. The last section, Section 6, concludes the paper.

Throughout the paper, we will use the following notations:(i) stands for the vector transpose.(ii)All vectors are interpreted as column vectors.(iii)Vectors are denoted by bold small letters.(iv) denotes the th element of vector .(v) represents the square of the Frobenius norm of .(vi) denotes a Gaussian distributed random variable with mean and variance .(vii) stands for expectation or sample mean of .(viii) stands for all integer numbers starting from to , thus .(ix) denotes the probability density function of .(x) is the speed of light.(xi) denotes the estimation of .(xii) stands for proportional.(xiii) defines the set for with .(xiv) denotes the uniform distribution on the interval .

2. Concept of Virtual Transmitters

Mathematically, the behavior of the multipath channel can be described by the time variant CIR , where indicates the discrete time instants and the delay [33]. According to [33], the CIR can be assumed to be constant for a short time interval at discrete time with index ,for , where is the number of MPCs, is the delay, the complex amplitude of the th MPC, and stands for the Dirac distribution [34] (please note that the CIR is generally a summation of an infinite number of MPCs; however, a practical receiver is only capable of capturing signals whose powers are above a certain sensitivity level). For notational conveniences, the LoS propagation path is considered also as a MPC in this paper. Assuming that the transmitted signal is band-limited with bandwidth and time-limited with a length smaller than , the signal received at time sampled with rate , bin indices , and the delay can be expressed aswhere denotes the white circular symmetric normal distributed receiver noise with variance . Using vector notation we obtain from (2)

In order to obtain the sparse structure of the CIR from the measurements , super resolution multipath estimation algorithms are necessary. The received signal is geometrically dependent on the transmitter and receiver positions as well as on the environment. Thus, the channel is spatially correlated as long as the spatial sampling is small enough. Hence, we use in this paper the dynamic multipath estimator named Kalman enhanced super resolution tracking (KEST) [20, 35–37] for estimating and tracking multipath parameters. KEST allows estimating the evolution of the CIR over time which is essential for Channel-SLAM as shown in the following section. KEST consists of a Kalman filter (KF) to estimate the complex amplitude and delay for each MPC utilizing maximum likelihood (ML) estimates as measurements. In the used implementation, KEST uses a standard model for the CIR which comprises a sum of weighted Dirac impulses as in (1). This model describes distinct paths sufficiently well. However, dense multipath components (DMC) lead to a model mismatch in the used KEST implementation. This model mismatch results in an increased variance of the estimated MPC parameters used as measurement noise in Channel-SLAM. For further details about KEST, see [20, 35–37].

To use the delay measurements of the tracked MPCs for positioning, a model describing the delays depending on the current user position is necessary. For developing such a model, we consider a static environment with a fixed transmitter and a receiver moving along an arbitrary trajectory. Figure 1 summarizes four propagation scenarios; for a detailed description see [31]. In the first scenario, the transmitted signal is reflected on a reflecting surface indicated by the blue lines. For reflection, we consider the effect of an electromagnetic wave reflected by a reflecting surface. When the receiver is moving, the reflection point on the reflecting surface is moving as well. If we mirror the physical transmitter position on the reflecting surface, we obtain the position of which is static during the receiver movement. The distance between and the receiver is equivalent to the propagation time of the reflected signal multiplied by the speed of light. Hence, the reflected signal can be interpreted as a direct signal from to the receiver.

Figure 1 

The figure shows four propagation scenarios: First scenario (blue): the transmitted signal is reflected on a reflecting surface. is defined by mirroring the physical transmitter position at the surface. Second scenario (red): the transmitted signal is reflected twice. is defined by mirroring the physical transmitter position at both surfaces. Third scenario (orange): the transmitted signal is scattered at . is defined at the position of . Fourth scenario (green): the transmitted signal is scattered and afterwards reflected on a reflecting surface. is defined by mirroring the scatterer at the surface.

This behavior can be extended to a multiple reflection scenario represented by the red lines. The transmitted signal is reflected two times. Equivalently, the location of can be determined by mirroring the transmitter position at both reflecting surfaces, as indicated in Figure 1. The distance between and the receiver is equivalent to the propagation time of the reflected signal multiplied by the speed of light. Thus, the signal reflected twice can also be interpreted as a direct signal from to the receiver.

Figure 1 exploits by the orange lines additionally a scenario where the signal is scattered, for example, at a lamp post. The propagation effect of scattering occurs if an electromagnetic wave impinges on an object and the energy is spread out in all directions [38]. Geometrically, the effect of scattering can be described as a fixed point at position in the pathway of the MPC. We define as at the position which is constant for all receiver positions for the MPC. Additionally, we treat , the constant distance between physical transmitter and scatterer, as an additional propagation distance associated with the MPC. Hence, the scattered signal can be interpreted as a direct signal from to the receiver, however, with a constant offset . Scattering and diffraction can be geometrically described as a fixed point at position in the pathway of the MPC and are considered as one model. Hence, unless otherwise stated, the description of scattering is equivalent for diffraction.

The fourth scenario considers the combination of both effects indicated in green. The transmitted signal is scattered at and afterwards reflected on the first reflecting surface. When the receiver is moving, the reflection point on the reflecting surface is moving as well. Hence, is defined by mirroring the scatterer at the first reflecting surface. Furthermore, between the transmitter and additional interactions are possible leading to the same position of .

To summarize, the propagation path of the th MPC can be equivalently described as a direct path with propagation length between and the receiver plus an additional constant propagation length ; hence,where denotes the speed of light and the position of the th VT (please note that the position of the VTs and the additional propagation lengths are constant over time. Nevertheless for notational convenience a time dependence on is introduced here). The additional propagation length is zero, that is, , if only reflections occurred on the pathway between physical transmitter and receiver or greater than zero, that is, , if the MPC is interacting with at least one scatterer. In general, can be interpreted as a clock offset between the th VT and the physical transmitter.

3. Channel-SLAM

3.1. Position Estimation

Figure 2 presents the available sensors together with the corresponding measurements. As shown on the left, we measure the sampled received signal as stated in (3) where we assume that the transmitter continuously emits known wideband signals. Based on , the multipath parameters amplitude and delay for each MPC are estimated and tracked by KEST. The estimated propagation path lengths of all MPCs of KEST are used as measurementsin Channel-SLAM with the corresponding variances . Because the VT positions are unknown, the receiver position and the positions of the VTs have to be estimated simultaneously. Thus, the state vector describing the complete system at time instant for MPCs iswith the receiver states and the VT states . The receiver state includes the receiver position , the receiver velocity , and the receiver’s clock bias ; hence,According to the description given in the previous section and (4), an MPC can be represented by a direct path between a VT and the receiver plus an additional propagation length. Hence, the parameters representing the th VT are defined aswhere is the position of the th VT and the additional propagation length. Using vector notation for all VTs, we obtain

Figure 2 

System model consisting of a terrestrial receiver and an IMU.

Additionally, as illustrated in Figure 2, an IMU is used. The IMU provides measurements of the acceleration and turn rates in three dimensions. After calibration, the heading change is used in Channel-SLAM as a control input and is therefore directly integrated into the transition model.

We use a discrete time representation for the transition and measurement model of the dynamic system withThe transition model in (10) describes the state evolution from time instant to time instant employing a possible nonlinear function with the process noise and using a control input which is in our case the heading change . The control input is considered as perfectly known and hence error-free. The measurement model (11) relates the state vector to the measurements by a possible nonlinear function and the measurement noise at time instant . Figure 3 shows the considered dynamic Bayesian network, that is, a first-order hidden Markov model.

Figure 3 

First-order hidden Markov model representing the dynamic system of Channel-SLAM.

Equations (10) and (11) can also be interpreted from a Bayesian perspective: based on measurements, we want to recursively estimate the unknown probability density function (PDF) of the state . In a recursive Bayesian formulation, this problem can be described as finding the posterior probability distributionRecursive Bayesian filtering provides a methodology to optimally estimate (12) by a prediction step to calculate and an update step to obtain which considers the measurement at time instant with the likelihood function [39, 40]. By assuming independence between the transition priors of the receiver state vector and the VT state vectors associated with the MPCs , the transition prior is defined here aswhere we inherently assume independence among MPCs, that is, propagation paths interacting with distinct objects. This is based on the well-known uncorrelated scattering assumption in wireless propagation channel modelling [38]. We obtain for the transition prior of the th MPCFor the transition prior of the receiver state vector we provide in Section 3.2 two models indicated by the function in Figure 3.

Assuming the elements of to be independent Gaussian distributed conditioned on the current state , can be expressed aswith the propagation lengthfor the th MPC, where denotes the corresponding variances.

3.2. Prediction Model Using Heading Changes

This paper considers a moving pedestrian carrying a hand-held device equipped with a terrestrial receiver and an IMU. A variety of pedestrian transition models exist in literature, for example, [41–44]; however, they do not fit for the considered application. Many of them focus on movements of groups, use additional information like floor plans, or do not incorporate information from an IMU. IMUs include in general accelerometers measuring acceleration and gyroscopes measuring turn rates , as indicated in Figure 2. These measurements are provided with respect to the sensor alignment [32], that is, the body frame. In order to obtain the measurements in a two-dimensional Cartesian coordinate system as shown in Figure 4, a transformation between the coordinate systems is necessary; see, for example, [45]. In our considered measurement scenario, the position of the IMU is assumed as constant with respect to the receiving antenna. Therefore, we are able to calculate the coordinate transformation matrices during a calibration phase when the pedestrian is standing still at the beginning. For other systems, where the sensor is decoupled, the sensor orientation has to be estimated continuously by applying strapdown navigation together with in-field calibration [46].

Figure 4 

Illustration of the prediction model for the pedestrian.

We propose two different constant velocity models, a linear model with Gaussian noise and a nonlinear model with Rician noise.

3.2.1. Gaussian-Transition-Model

The first proposed transition model is based on a discrete white noise acceleration model [47], referred to as Gaussian-Transition-Model, within a two-dimensional Cartesian coordinate system. The receiver state vector consists of the - positionsand the velocitieswhere are the corresponding velocities in - direction and the receiver’s clock bias where a standard clock bias model is used [12, 48] (Please note that if transmitter and receiver oscillators provide different frequencies, a clock drift parameter has to be considered additionally). The transition matrix in (17) includes a rotation matrix with the heading changes , with where and is the transition noise of the receiver state vector with covariancewhere defines the continuous-time process noise intensity that has to be set based on the application with physical dimension and the variance of the clock bias. This transition model is similar to the transition model presented in [29, 30], with the advantage that the transition model is linear if the heading changes are known.

Since we incorporate only the heading changes , we do not have any speed measurements and the speed has to be estimated implicitly. In order to adapt quickly to different walking speeds, has to be large or many particles have to be used to cover all possible movements. A large value for may cause backward movements of the transition model which results in estimation errors of Channel-SLAM. Another drawback of this transition model is that the estimated state information is completely lost when the user is standing; thus, the velocity components are zero. In order to overcome the mentioned problems, we develop a second transition model as described in the following.

3.2.2. Rician-Transition-Model

Similarly to Section 3.2.1, we follow a two-dimensional positioning approach in the Cartesian coordinate system with the receiver state vector which consists of the - receiver positions as defined in (18), the receiver velocity vector , and the receiver’s clock bias . The receiver velocity vector iswhere is the receiver speed and the heading of the receiver; see Figure 4. The heading describes the walking direction of the pedestrian with respect to the Cartesian coordinate system. Hence, we can define the transition model withwhere the velocity follows a Rician distribution withwith scale parameter . For speeds close to zero, the Rician distribution approximates a Rayleigh distribution; thus, the speed is always positive. Hence, it has the advantage of preventing the filter from converging to negative velocities, which are highly unlikely regarding a normal pedestrian walking behavior (Please note that the developed transition model does not include standing or walking backwards phases. This could be additionally considered by extracting more information from the IMU measurements). This is important for our approach, since ambiguities of Channel-SLAM could otherwise cause a movement in the wrong direction. On the other hand, for higher speeds, the distribution becomes approximately Gaussian which reflects empirical data of pedestrian walking speeds [43].

Finally, the heading of the user is defined bywhere is the heading change from the IMU after calibration with the heading noise using a von Mises distribution. For the transition prior of the clock bias we use similar to Section 3.2.1 a standard clock bias model , where defines the transition noise with the variance [12, 48] (Please note that if transmitter and receiver oscillators provide different frequencies, a clock drift parameter has to be considered additionally). In the following we refer to this transition model as Rician-Transition-Model.

3.3. Rao-Blackwellized Particle Filter

As introduced in [31], Channel-SLAM is derived based on Rao-Blackwellization where the state space of is partitioned into subspaces. Hence, we use particle filters (PFs) to estimate the subspaces representing the VTs inside a PF. The reason to use a PF instead of a low complexity EKF is the high nonlinearity of the measurements in (16). As shown in Figure 5, the algorithm is based on a superordinate particle filter (superPF) and subordinate particle filters (subPFs): Each particle of the superPF with the state vector consists of subPFs. Each subPF is represented by the particles with , where stands for the number of particles in the th subPF with , estimating . Using subPFs for each VT allows using different numbers of particles in each subPF and, furthermore, allows a dynamic adjustment of the number of particles for each subPF introduced in Section 3.4.

Figure 5 

The algorithm is based on a superordinate particle filter (superPF) and subordinate particle filters (subPFs). Each particle of the superPF consists of subPFs.

According to [31], the marginalized posterior filtered density of the superPF can be approximated by importance samples (see [31, 39]) aswhere defines the weight for the th particle at time instant withand the weight of the subPFs at time instant withContrarily to [31], resampling is performed at each time instant to prevent degeneration; hence, (27) and (28) do not depend on the weights and . Additionally, the derivations in [31] consider a regularized PF [39] where is drawn after resampling from the Gaussian-Kernel . The Gaussian-Kernel improves the robustness of Channel-SLAM to cope with small model mismatches in the measurements.

3.4. Particle Filter Implementation

Algorithm 1 provides the pseudocode of Channel-SLAM, which is executed at every time instant with the estimates obtained from KEST. During the initialization, at time instant , the particles of the superPF are initialized according to prior knowledge. The particles of the subPFs are initialized dependent on and the measurements for the th MPC. To initialize the states of with of the th subPF associated with the th MPC a grid is used. The positions of the particles are distributed on a grid inside a circle around with radius such thatwith spacing (the number of grid points can be estimated by Gauss’s circle problem). The additional propagation length is , where we inherently assume for the initialization. Hence, the total number of particles can be calculated asFor each particle of the superPF, the receiver state is propagated according to the transition model described in Section 3.2 indicated by the function using the heading changes (cf. Line () in Algorithm 1). Afterwards, Channel-SLAM determines whether the number of tracked MPCs has changed. In case that new MPCs have been detected, new subPFs are added and initialized using (29) (cf. Line () in Algorithm 1). In case that MPCs are not tracked by KEST anymore, the corresponding subPFs are removed (cf. Line () in Algorithm 1). Equivalent to [31], neither KEST nor Channel-SLAM considers retracking of previous MPCs. Hence, if the tracking of an MPC has been lost and is regained again, the corresponding VT is initialized without any prior information. According to (14), the state is time-invariant; hence, each subPF assigns the states of the VTs with . Thereafter, the weights for the subPFs and superPF are calculated using (28) and (27).