# Positioning Using Terrestrial Multipath Signals and Inertial Sensors

Authors: Christian Gentner; Robert Pöhlmann; Markus Ulmschneider; Thomas Jost; Siwei Zhang
Journal: Mobile Information Systems
 (2017)

---

## 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.


---

## Body

## 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: Section2 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)
[·]T stands for the vector transpose.(ii)
All vectors are interpreted as column vectors.(iii)
Vectors are denoted by bold small letters.(iv)
xl denotes the lth element of vector x.(v)
A2=∑l∑mAl,m2 represents the square of the Frobenius norm of A.(vi)
a~Nμa,σa2 denotes a Gaussian distributed random variable a with mean μa and variance σa2.(vii)
Ex stands for expectation or sample mean of x.(viii)
1:k stands for all integer numbers starting from 1 to k, thus 1,2,…,k.(ix)
p(x) denotes the probability density function of x.(x)
c is the speed of light.(xi)
x^ denotes the estimation of x.(xii)
∝ stands for proportional.(xiii)
{x(i)}i=1N defines the set for xi with i=1,…,N.(xiv)
U0,N denotes the uniform distribution on the interval [0,N].


## 2. Concept of Virtual Transmitters

Mathematically, the behavior of the multipath channel can be described by the time variant CIRh(tk,τ), where tk indicates the discrete time instants and τ the delay [33]. According to [33], the CIR h(tk,τ) can be assumed to be constant for a short time interval T at discrete time tk with index k,(1)htk,τ=∑i=0Ntk-1αitk·δτ-τitk,for T0≤tk≤T0+T, where N(tk) is the number of MPCs, τi(tk) is the delay, αi(tk) the complex amplitude of the ith 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 s(tk) is band-limited with bandwidth B and time-limited with a length smaller than T, the signal received at time tk sampled with rate B, bin indices m=0,…,M-1, and the delay τm=m/B can be expressed as(2)ytk,τm=∑i=0Ntk-1αitksτm-τitk+nτm,where n(τm) denotes the white circular symmetric normal distributed receiver noise with variance σn2. Using vector notation we obtain from (2)(3)ytk=ytk,τ0,…,ytk,τm,…,ytk,τM-1.In order to obtain the sparse structure of the CIR from the measurementsy(tk), 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 α^i(tk) and delay τ^i(tk) for each MPC i 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τi(tk) depending on the current user position ru(tk) 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 rVT,1 of VT1 which is static during the receiver movement. The distance between VT1 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 VT1 to the receiver.Figure 1
The figure shows four propagation scenarios: First scenario (blue): the transmitted signal is reflected on a reflecting surface.VT1 is defined by mirroring the physical transmitter position at the surface. Second scenario (red): the transmitted signal is reflected twice. VT2 is defined by mirroring the physical transmitter position at both surfaces. Third scenario (orange): the transmitted signal is scattered at S. VT3 is defined at the position of S. Fourth scenario (green): the transmitted signal is scattered and afterwards reflected on a reflecting surface. VT4 is defined by mirroring the scatterer S 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 ofVT2 can be determined by mirroring the transmitter position at both reflecting surfaces, as indicated in Figure 1. The distance between VT2 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 VT2 to the receiver.Figure1 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 S at position rS in the pathway of the MPC. We define S as VT3 at the position rS which is constant for all receiver positions for the MPC. Additionally, we treat dVT>0, 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 VT3 to the receiver, however, with a constant offset dVT. Scattering and diffraction can be geometrically described as a fixed point S at position rS 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 atS 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, VT4 is defined by mirroring the scatterer S at the first reflecting surface. Furthermore, between the transmitter and S additional interactions are possible leading to the same position of VT4.To summarize, the propagation path of theith MPC can be equivalently described as a direct path with propagation length di(tk) between VTi and the receiver plus an additional constant propagation length dVT,i(tk); hence,(4)ditk=τitk·c=rutk-rVT,itk+dVT,itk,where c denotes the speed of light and rVT,i(tk) the position of the ith 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 tk is introduced here). The additional propagation length is zero, that is, dVT,i(tk)=0, if only reflections occurred on the pathway between physical transmitter and receiver or greater than zero, that is, dVT,i(tk)>0, if the MPC is interacting with at least one scatterer. In general, dVT,i(tk)/c can be interpreted as a clock offset between the ith VT and the physical transmitter.


## 3. Channel-SLAM

### 3.1. Position Estimation

Figure2 presents the available sensors together with the corresponding measurements. As shown on the left, we measure the sampled received signal y(tk) as stated in (3) where we assume that the transmitter continuously emits known wideband signals. Based on y(tk), the multipath parameters amplitude αi(tk) and delay τi(tk)=di(tk)/c for each MPC are estimated and tracked by KEST. The estimated propagation path lengths d^i(tk)=τ^i(tk)·c of all N(tk) MPCs of KEST are used as measurements(5)ztk=d^0tk,…,d^Ntk-1tkTin Channel-SLAM with the corresponding variances σz(tk). Because the VT positions are unknown, the receiver position and the positions of the VTs have to be estimated simultaneously. Thus, the state vector x(tk) describing the complete system at time instant tk for N(tk) MPCs is(6)xtk=xutkT,xVTtkTT,with the receiver states xu(tk) and the VT states xVT(tk). The receiver state xu(tk) includes the receiver position ru(tk), the receiver velocity vu(tk), and the receiver’s clock bias bu(tk); hence,(7)xutk=rutkT,vutkT,butkT.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 ith VT are defined as(8)xVT,itk=rVT,itkT,dVT,itkT,where rVT,i(tk) is the position of the ith VT and dVT,i(tk) the additional propagation length. Using vector notation for all VTs, we obtain(9)xVTtk=xVT,0tkT,…,xVT,Ntk-1tkTT.Figure 2
System model consisting of a terrestrial receiver and an IMU.Additionally, as illustrated in Figure2, an IMU is used. The IMU provides measurements of the acceleration ab(tk) and turn rates ωibb(tk) in three dimensions. After calibration, the heading change Ψ˙(tk) 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 with(10)xtk=fxtk-1,Ψ˙tk,nttk,(11)ztk=hxtk,nhtk.The transition model in (10) describes the state evolution from time instant tk-1 to time instant tk employing a possible nonlinear function f(·,·,·) with the process noise nt(tk) and using a control input which is in our case the heading change Ψ˙(tk). 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 h(·,·) and the measurement noise nh(tk) at time instant tk. 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 x(tk). In a recursive Bayesian formulation, this problem can be described as finding the posterior probability distribution(12)pxtk∣zt1:k,Ψ˙t1:k,xt0.Recursive Bayesian filtering provides a methodology to optimally estimate (12) by a prediction step to calculate pxtk∣zt1:k-1,Ψ˙t1:k,xt0 and an update step to obtain p(xtk∣zt1:k,Ψ˙t1:k,xt0) which considers the measurement z(tk) at time instant tk with the likelihood function p(z(tk)∣xtk) [39, 40]. By assuming independence between the transition priors of the receiver state vector xu(tk) and the VT state vectors xVT,i(tk) associated with the MPCs i=0,…,N(tk)-1, the transition prior pxtk∣xtk-1,Ψ˙tk,xt0 is defined here as(13)pxtk∣xtk-1,Ψ˙tk,xt0=pxutk∣xutk-1,Ψ˙tk,xut0×∏i=0Ntk-1pxVT,itk∣xVT,itk-1,where 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 pxVT,i(tk)∣xVT,itk-1 of the ith MPC(14)pxVT,itk∣xVT,itk-1=δxVT,itk-xVT,itk-1.For the transition prior p(xu(tk)∣xutk-1,Ψ˙tk) of the receiver state vector we provide in Section 3.2 two models indicated by the function f(xu(tk-1),Ψ˙(tk),nu(tk)) in Figure 3.Assuming the elements ofz(tk) to be independent Gaussian distributed conditioned on the current state x(tk), p(z(tk)∣xtk) can be expressed as(15)pztk∣xtk=∏i=0Ntk-112πσd,itke-d^itk-ditk2/2σd,i2tkwith the propagation length(16)ditk=rutk-rVT,itk+dVT,itk+butk·c,for the ith MPC, where σd,i2(tk) 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 ab(tk) and gyroscopes measuring turn rates ωibb(tk), 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, with(17)x-utk=Autδ,Ψ˙tkx-utk-1+nttk,in a two-dimensional Cartesian coordinate system. The receiver state vector x-u(tk)=[r-u(tk),v-u(tk),bu(tk)]T consists of the x-y positions(18)r-utk=ru,xtk,ru,ytkTand the velocities(19)v-utk=vu,xtk,vu,ytkT,where vu,x(tk),vu,y(tk) are the corresponding velocities in x-y direction and the receiver’s clock bias bu(tk) 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 Au(tδ,Ψ˙(tk)) in (17) includes a rotation matrix with the heading changes Ψ˙(tk), with (20)Autδ,Ψ˙tk=10tδ00010tδ000cos⁡Ψ˙tk-sin⁡Ψ˙tk000sin⁡Ψ˙tkcos⁡Ψ˙tk000001,where tδ=tk-tk-1 and nt(tk)~N0,Qu(tδ) is the transition noise of the receiver state vector with covariance(21)Qutδ=σqu2tδ330σqu2tδ22000σqu2tδ330σqu2tδ220σqu2tδ220σqu2tδ000σqu2tδ220σqu2tδ00000σqb2,where σqu2 defines the continuous-time process noise intensity that has to be set based on the application with physical dimension [m2/s3] and σqb2 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Ψ˙(tk), we do not have any speed measurements and the speed has to be estimated implicitly. In order to adapt quickly to different walking speeds, σqu has to be large or many particles have to be used to cover all possible movements. A large value for σqu 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 Section3.2.1, we follow a two-dimensional positioning approach in the Cartesian coordinate system with the receiver state vector x-u(tk)=[r˘u(tk),v˘u(tk),bu(tk)]T which consists of the x-y receiver positions r˘u(tk) as defined in (18), the receiver velocity vector v˘u(tk), and the receiver’s clock bias bu(tk). The receiver velocity vector is(22)v˘utk=vutk,ΨutkT,where vu(tk) is the receiver speed and Ψu(tk) the heading of the receiver; see Figure 4. The heading Ψu(tk) describes the walking direction of the pedestrian with respect to the Cartesian coordinate system. Hence, we can define the transition model with(23)r˘utk=r˘utk-1+tδvutkcos⁡Ψutksin⁡Ψutk,where the velocity follows a Rician distribution with(24)vutk~Rvutk-1,σmtk,with scale parameter σm(tk). 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 by(25)Ψutk=Ψutk-1+tk-tk-1Ψ˙tk+wΨtk,where Ψ˙(tk) is the heading change from the IMU after calibration with the heading noise wΨ(tk) 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 bu(tk)=bu(tk-1)+nb(tk-1), where nb(tk-1) defines the transition noise with the variance σqb2 [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 x(tk) 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 j=1⋯Ns of the superPF with the state vector xu(j)(tk)=ru(j)(tk)T,vu(j)(tk)T,bu(j)(tk)T consists of N(tk) subPFs. Each subPF is represented by the particles xVT,i(j,a)(tk) with a=1,…,NP,j,i(tk), where NP,j,i(tk) stands for the number of particles in the ith subPF with i=0,…,N(tk)-1, estimating xVT,i(j)(tk). 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 particlej=1⋯Ns of the superPF consists of N(tk) subPFs.According to [31], the marginalized posterior filtered density p(xutk∣zt1:k,Ψ˙t1:k) of the superPF can be approximated by importance samples (see [31, 39]) as(26)pxutk∣zt1:k,Ψ˙t1:k≈∑j=1Nswjtkδxutk-xujtk,where w(j)(tk) defines the weight for the jth particle at time instant tk with(27)wjtk∝pztk∣xujtk,ztk-1∝∏i=0Nt-1∑a=1NP,i,jtkwij,atkand the weight wi(j,a)(tk) of the subPFs at time instant tk with(28)wij,atk≜pd^itk∣xujtk,xVT,ij,atk.Contrarily to [31], resampling is performed at each time instant to prevent degeneration; hence, (27) and (28) do not depend on the weights w(j)(tk-1) and wi(j,a)(tk-1). Additionally, the derivations in [31] consider a regularized PF [39] where xVT,i(j,a)(tk) is drawn after resampling from the Gaussian-Kernel K(·). The Gaussian-Kernel K(·) improves the robustness of Channel-SLAM to cope with small model mismatches in the measurements.


### 3.4. Particle Filter Implementation

Algorithm1 provides the pseudocode of Channel-SLAM, which is executed at every time instant tk≥t0 with the estimates z(tk),σz(tk) obtained from KEST. During the initialization, at time instant tk=t0, the particles {xu(j)(t0)}j=1Ns of the superPF are initialized according to prior knowledge. The particles {xVT,i(j,a)(t0)}a=1Np,j,i of the subPFs are initialized dependent on xu(j)(t0) and the measurements d^i(t0) for the ith MPC. To initialize the states of xVT,i(j,a)(t0) with a=1,…,Np,j,i of the jth subPF associated with the ith MPC a grid is used. The positions rVT,i(j,a)(t0) of the particles {xVT,i(j,a)(t0)}a=1NP,j,i are distributed on a grid inside a circle around ru(j)(t0) with radius d^i(t0)+Δd such that(29)0≤rVT,ij,at0-rujt0≤d^it0+Δdwith spacing Δd (the number of grid points NP,j,i(tk) can be estimated by Gauss’s circle problem). The additional propagation length is dVT,i(j,a)(t0)=d^i(t0)-rVT,i(j,a)(t0)-ru(j)(t0), where we inherently assume bu(t0)=0 for the initialization. Hence, the total number of particles can be calculated as(30)Nttk=∑j=1Np∑i=0Ntk-1NP,j,itk.For each particle j of the superPF, the receiver state xu(j)(tk-1) is propagated according to the transition model described in Section 3.2 indicated by the function f(xu(j)(tk-1),Ψ˙(tk),nt(tk)) using the heading changes Ψ˙(tk) (cf. Line (5) 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 (7) in Algorithm 1). In case that MPCs are not tracked by KEST anymore, the corresponding subPFs are removed (cf. Line (9) 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 xVT,i(tk) is time-invariant; hence, each subPF assigns the states of the VTs with xVT,i(j,a)(tk)=xVT,i(j,a)(tk-1). Thereafter, the weights for the subPFs and superPF are calculated using (28) and (27).Algorithm 1:Channel-SLAM for time instant tk.
Input:
Multipath estimates: z(tk), σz(tk) and Ψ˙(tk);
States for tk>t0:
x
u
(
j
)
(
t
k
-
1
)
,
x
V
T
,
i
(
j
,
a
)
(
t
k
-
1
)
a
=
1
N
p
,
j
,
i
i
=
0
N
(
t
k
-
1
)
-
1
j
=
1
N
s
Output:
States for ≥t0:
x
u
(
j
)
(
t
k
)
,
x
V
T
,
i
(
j
,
a
)
(
t
k
)
a
=
1
N
p
,
j
,
i
i
=
0
N
(
t
k
)
-
1
j
=
1
N
s
MMSE estimate: x^(tk) for tk>t0
(1) iftk=t0then
(2) Initialization using z(tk) and σz(tk);
(3) else
(4) forj=1:Nsdo
(5) Draw xu(j)(tk)=f(xu(j)(tk-1),Ψ˙(tk),w(tk-1));
(6) ifNew paths detectedthen
(7) Initialize new subPFs;
(8) ifTracking of paths lostthen
(9) Delete corresponding subPFs;
(10) fori=0:N(tk)-1do
(11) fora=1:NP,j,i(tk)do
(12) Assign xVT,i(j,a)(tk)=xVT,i(j,a)(tk-1);
(13) Calculate wij,a(tk)=p(zi(tk)∣xujtk,xVT,ij,atk);
(14) Calculate total subPF weight:
t
j
,
i
=
S
U
M
[
{
w
i
(
j
,
a
)
(
t
k
)
}
a
=
1
N
P
,
j
,
i
(
t
k
)
];
(15) w(j)(tk)=∏i=0N(tk)-1‍tj,i;
(16) Resample using Algorithm 2;
(17) Calculate MMSE x^(tk) according to (31);Algorithm 2:Channel-SLAM resampling algorithm.
Input:
States and weights:
x
u
j
(
t
k
)
,
w
(
j
)
(
t
k
)
,
x
V
T
,
i
(
j
,
a
)
(
t
k
)
,
w
i
(
j
,
a
)
(
t
k
)
a
=
1
N
p
,
j
,
i
(
t
k
)
i
=
0
N
(
t
k
)
-
1
j
=
1
N
s
Output:
Resampled states and weights:
x
~
u
j
(
t
k
)
,
w
~
j
(
t
k
)
,
x
~
V
T
,
i
j
,
a
(
t
k
)
,
w
~
i
(
j
,
a
)
(
t
k
)
a
=
1
N
p
,
j
,
i
(
t
k
)
i
=
0
N
(
t
k
)
-
1
j
=
1
N
s
(1) Initialize the CDF for superPF: cp1=w1(tk);
(2) forj=2:Nsdo
(3) Construct CDF for superPF:
c
p
j
=
c
p
j
-
1
+
w
j
(
t
k
);
(4) f=1;
(5) Draw starting point: up1~U0,Ns-1;
(6) forj=1:Nsdo
(7) upj=up1+Ns-1(j-1);
(8) whileupj>cpfdo
(9) f=f+1;
(10) Assign: x~uj(tk),w~j(tk)=xuftk,1/Ns;
(11) fori=0:N(tk)-1do
(12) Initialize the CDF for (f,i)-th subPF:
c
b
1
=
w
i
(
f
,
1
)
(
t
k
);
(13) fora=2:NP,i,f(tk)do
(14) Construct CDF for subPF:
c
b
a
=
c
b
a
-
1
+
w
i
(
f
,
a
)
(
t
k
);
(15) g=1,ar=1;
(16) fora=1:NP,f,i(tk)do
(17) uba=ub1+1NP,i,f(tk)(a-1);
(18) whileuba>cbgdo
(19) g=g+1;
(20) ifB(x~VT,ij,artk)≤Nmthen
(21) Assign: x~VT,ij,ar(tk)=xVT,if,g(tk);
(22) ar=ar+1;
(23) Bx~VT,ij,artk=Bx~VT,ij,artk+1;
(24) Update number of particles: NP,j,i(tk)=ar;
(25) fora=1:NP,j,i(tk)do
(26) Draw x~VT,ij,atk from the Gaussian-Kernel;
(27) Assign: w~ij,a(tk)=1/NP,j,i(tk);Afterwards, the subPFs and superPF are resampled. The basic idea of the resampling method is to eliminate particles with low weights and reproduce particles with high weights. Algorithm2 shows a pseudocode of the resampling algorithm of Channel-SLAM which is based on the systematic resampling algorithm [49]. Similarly to Algorithm 1, the Channel-SLAM resampling algorithm consists of a resampling algorithm for the superPF which includes resampling algorithms for the subPFs. First, the estimated sampled cumulative distribution function (CDF) of the superPF is constructed, presented by a vector cp with length Np and element cpj with j=1,…,Np. According to the estimated sampled CDF of the superPF, the subPFs are eliminated or resampled. The particles of the superPF with index f are assigned to the resampled particle index j; see Algorithm 2, Line (10), for the assignment of the receiver state. Afterwards, the (f,i)th subPF is resampled with i=0,…,N(tk)-1 using a systematic resampling algorithm, where cb represents the estimated sampled CDF of the (f,i)th subPF.As mentioned before, whenever a new MPC is tracked, many particles are initialized to cover all possible VT positions in each subPF. For example if theith MPC has a delay of d^i(t0) = 30 m, each ith subPF would use NP,j,i(t0)=2821 particles according to (29) for j=1,…,Np with spacing Δd = 1 m. However during the receiver movement many particles of the subPFs are resampled at the same grid point because the states of the VTs xVT,i(tk) are time-invariant. In order to adapt the number of particles, we limit the number of resampled particles per grid point to Nm; see Algorithm 2, Line (20). Evaluations in Section 5 show that the reduction of the number of particles does not influence the positioning accuracy but, however, leads to a gain on computational performance. Afterwards, the states of the VTs xVT,i(j,a)(tk) are drawn using a Gaussian-Kernel (cf. Line (26) in Algorithm 2). As stated in [31], the Gaussian-Kernel has a low bandwidth which does not influence the grid based reduction method.Usually, we are interested in a point estimate of the state instead of its a posteriori PDF. According to [31], the MMSE point estimate, x^(tk)=x^u(tk)T,x^VT(tk)TT (cf. Line (17) in Algorithm 1), is calculated by(31)x^tk≈∑j=1Nswjtk×xujtk∑a=1NP,j,0tkw0j,atkxVT,0j,atk⋮∑a=1NP,j,Ntk-1tkwNtk-1j,atkxVT,Ntk-1j,atk.


## 3.1. Position Estimation

Figure2 presents the available sensors together with the corresponding measurements. As shown on the left, we measure the sampled received signal y(tk) as stated in (3) where we assume that the transmitter continuously emits known wideband signals. Based on y(tk), the multipath parameters amplitude αi(tk) and delay τi(tk)=di(tk)/c for each MPC are estimated and tracked by KEST. The estimated propagation path lengths d^i(tk)=τ^i(tk)·c of all N(tk) MPCs of KEST are used as measurements(5)ztk=d^0tk,…,d^Ntk-1tkTin Channel-SLAM with the corresponding variances σz(tk). Because the VT positions are unknown, the receiver position and the positions of the VTs have to be estimated simultaneously. Thus, the state vector x(tk) describing the complete system at time instant tk for N(tk) MPCs is(6)xtk=xutkT,xVTtkTT,with the receiver states xu(tk) and the VT states xVT(tk). The receiver state xu(tk) includes the receiver position ru(tk), the receiver velocity vu(tk), and the receiver’s clock bias bu(tk); hence,(7)xutk=rutkT,vutkT,butkT.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 ith VT are defined as(8)xVT,itk=rVT,itkT,dVT,itkT,where rVT,i(tk) is the position of the ith VT and dVT,i(tk) the additional propagation length. Using vector notation for all VTs, we obtain(9)xVTtk=xVT,0tkT,…,xVT,Ntk-1tkTT.Figure 2
System model consisting of a terrestrial receiver and an IMU.Additionally, as illustrated in Figure2, an IMU is used. The IMU provides measurements of the acceleration ab(tk) and turn rates ωibb(tk) in three dimensions. After calibration, the heading change Ψ˙(tk) 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 with(10)xtk=fxtk-1,Ψ˙tk,nttk,(11)ztk=hxtk,nhtk.The transition model in (10) describes the state evolution from time instant tk-1 to time instant tk employing a possible nonlinear function f(·,·,·) with the process noise nt(tk) and using a control input which is in our case the heading change Ψ˙(tk). 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 h(·,·) and the measurement noise nh(tk) at time instant tk. 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 x(tk). In a recursive Bayesian formulation, this problem can be described as finding the posterior probability distribution(12)pxtk∣zt1:k,Ψ˙t1:k,xt0.Recursive Bayesian filtering provides a methodology to optimally estimate (12) by a prediction step to calculate pxtk∣zt1:k-1,Ψ˙t1:k,xt0 and an update step to obtain p(xtk∣zt1:k,Ψ˙t1:k,xt0) which considers the measurement z(tk) at time instant tk with the likelihood function p(z(tk)∣xtk) [39, 40]. By assuming independence between the transition priors of the receiver state vector xu(tk) and the VT state vectors xVT,i(tk) associated with the MPCs i=0,…,N(tk)-1, the transition prior pxtk∣xtk-1,Ψ˙tk,xt0 is defined here as(13)pxtk∣xtk-1,Ψ˙tk,xt0=pxutk∣xutk-1,Ψ˙tk,xut0×∏i=0Ntk-1pxVT,itk∣xVT,itk-1,where 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 pxVT,i(tk)∣xVT,itk-1 of the ith MPC(14)pxVT,itk∣xVT,itk-1=δxVT,itk-xVT,itk-1.For the transition prior p(xu(tk)∣xutk-1,Ψ˙tk) of the receiver state vector we provide in Section 3.2 two models indicated by the function f(xu(tk-1),Ψ˙(tk),nu(tk)) in Figure 3.Assuming the elements ofz(tk) to be independent Gaussian distributed conditioned on the current state x(tk), p(z(tk)∣xtk) can be expressed as(15)pztk∣xtk=∏i=0Ntk-112πσd,itke-d^itk-ditk2/2σd,i2tkwith the propagation length(16)ditk=rutk-rVT,itk+dVT,itk+butk·c,for the ith MPC, where σd,i2(tk) 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 ab(tk) and gyroscopes measuring turn rates ωibb(tk), 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, with(17)x-utk=Autδ,Ψ˙tkx-utk-1+nttk,in a two-dimensional Cartesian coordinate system. The receiver state vector x-u(tk)=[r-u(tk),v-u(tk),bu(tk)]T consists of the x-y positions(18)r-utk=ru,xtk,ru,ytkTand the velocities(19)v-utk=vu,xtk,vu,ytkT,where vu,x(tk),vu,y(tk) are the corresponding velocities in x-y direction and the receiver’s clock bias bu(tk) 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 Au(tδ,Ψ˙(tk)) in (17) includes a rotation matrix with the heading changes Ψ˙(tk), with (20)Autδ,Ψ˙tk=10tδ00010tδ000cos⁡Ψ˙tk-sin⁡Ψ˙tk000sin⁡Ψ˙tkcos⁡Ψ˙tk000001,where tδ=tk-tk-1 and nt(tk)~N0,Qu(tδ) is the transition noise of the receiver state vector with covariance(21)Qutδ=σqu2tδ330σqu2tδ22000σqu2tδ330σqu2tδ220σqu2tδ220σqu2tδ000σqu2tδ220σqu2tδ00000σqb2,where σqu2 defines the continuous-time process noise intensity that has to be set based on the application with physical dimension [m2/s3] and σqb2 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Ψ˙(tk), we do not have any speed measurements and the speed has to be estimated implicitly. In order to adapt quickly to different walking speeds, σqu has to be large or many particles have to be used to cover all possible movements. A large value for σqu 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 Section3.2.1, we follow a two-dimensional positioning approach in the Cartesian coordinate system with the receiver state vector x-u(tk)=[r˘u(tk),v˘u(tk),bu(tk)]T which consists of the x-y receiver positions r˘u(tk) as defined in (18), the receiver velocity vector v˘u(tk), and the receiver’s clock bias bu(tk). The receiver velocity vector is(22)v˘utk=vutk,ΨutkT,where vu(tk) is the receiver speed and Ψu(tk) the heading of the receiver; see Figure 4. The heading Ψu(tk) describes the walking direction of the pedestrian with respect to the Cartesian coordinate system. Hence, we can define the transition model with(23)r˘utk=r˘utk-1+tδvutkcos⁡Ψutksin⁡Ψutk,where the velocity follows a Rician distribution with(24)vutk~Rvutk-1,σmtk,with scale parameter σm(tk). 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 by(25)Ψutk=Ψutk-1+tk-tk-1Ψ˙tk+wΨtk,where Ψ˙(tk) is the heading change from the IMU after calibration with the heading noise wΨ(tk) 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 bu(tk)=bu(tk-1)+nb(tk-1), where nb(tk-1) defines the transition noise with the variance σqb2 [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.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, with(17)x-utk=Autδ,Ψ˙tkx-utk-1+nttk,in a two-dimensional Cartesian coordinate system. The receiver state vector x-u(tk)=[r-u(tk),v-u(tk),bu(tk)]T consists of the x-y positions(18)r-utk=ru,xtk,ru,ytkTand the velocities(19)v-utk=vu,xtk,vu,ytkT,where vu,x(tk),vu,y(tk) are the corresponding velocities in x-y direction and the receiver’s clock bias bu(tk) 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 Au(tδ,Ψ˙(tk)) in (17) includes a rotation matrix with the heading changes Ψ˙(tk), with (20)Autδ,Ψ˙tk=10tδ00010tδ000cos⁡Ψ˙tk-sin⁡Ψ˙tk000sin⁡Ψ˙tkcos⁡Ψ˙tk000001,where tδ=tk-tk-1 and nt(tk)~N0,Qu(tδ) is the transition noise of the receiver state vector with covariance(21)Qutδ=σqu2tδ330σqu2tδ22000σqu2tδ330σqu2tδ220σqu2tδ220σqu2tδ000σqu2tδ220σqu2tδ00000σqb2,where σqu2 defines the continuous-time process noise intensity that has to be set based on the application with physical dimension [m2/s3] and σqb2 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Ψ˙(tk), we do not have any speed measurements and the speed has to be estimated implicitly. In order to adapt quickly to different walking speeds, σqu has to be large or many particles have to be used to cover all possible movements. A large value for σqu 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 Section3.2.1, we follow a two-dimensional positioning approach in the Cartesian coordinate system with the receiver state vector x-u(tk)=[r˘u(tk),v˘u(tk),bu(tk)]T which consists of the x-y receiver positions r˘u(tk) as defined in (18), the receiver velocity vector v˘u(tk), and the receiver’s clock bias bu(tk). The receiver velocity vector is(22)v˘utk=vutk,ΨutkT,where vu(tk) is the receiver speed and Ψu(tk) the heading of the receiver; see Figure 4. The heading Ψu(tk) describes the walking direction of the pedestrian with respect to the Cartesian coordinate system. Hence, we can define the transition model with(23)r˘utk=r˘utk-1+tδvutkcos⁡Ψutksin⁡Ψutk,where the velocity follows a Rician distribution with(24)vutk~Rvutk-1,σmtk,with scale parameter σm(tk). 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 by(25)Ψutk=Ψutk-1+tk-tk-1Ψ˙tk+wΨtk,where Ψ˙(tk) is the heading change from the IMU after calibration with the heading noise wΨ(tk) 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 bu(tk)=bu(tk-1)+nb(tk-1), where nb(tk-1) defines the transition noise with the variance σqb2 [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 x(tk) 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 j=1⋯Ns of the superPF with the state vector xu(j)(tk)=ru(j)(tk)T,vu(j)(tk)T,bu(j)(tk)T consists of N(tk) subPFs. Each subPF is represented by the particles xVT,i(j,a)(tk) with a=1,…,NP,j,i(tk), where NP,j,i(tk) stands for the number of particles in the ith subPF with i=0,…,N(tk)-1, estimating xVT,i(j)(tk). 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 particlej=1⋯Ns of the superPF consists of N(tk) subPFs.According to [31], the marginalized posterior filtered density p(xutk∣zt1:k,Ψ˙t1:k) of the superPF can be approximated by importance samples (see [31, 39]) as(26)pxutk∣zt1:k,Ψ˙t1:k≈∑j=1Nswjtkδxutk-xujtk,where w(j)(tk) defines the weight for the jth particle at time instant tk with(27)wjtk∝pztk∣xujtk,ztk-1∝∏i=0Nt-1∑a=1NP,i,jtkwij,atkand the weight wi(j,a)(tk) of the subPFs at time instant tk with(28)wij,atk≜pd^itk∣xujtk,xVT,ij,atk.Contrarily to [31], resampling is performed at each time instant to prevent degeneration; hence, (27) and (28) do not depend on the weights w(j)(tk-1) and wi(j,a)(tk-1). Additionally, the derivations in [31] consider a regularized PF [39] where xVT,i(j,a)(tk) is drawn after resampling from the Gaussian-Kernel K(·). The Gaussian-Kernel K(·) improves the robustness of Channel-SLAM to cope with small model mismatches in the measurements.


## 3.4. Particle Filter Implementation

Algorithm1 provides the pseudocode of Channel-SLAM, which is executed at every time instant tk≥t0 with the estimates z(tk),σz(tk) obtained from KEST. During the initialization, at time instant tk=t0, the particles {xu(j)(t0)}j=1Ns of the superPF are initialized according to prior knowledge. The particles {xVT,i(j,a)(t0)}a=1Np,j,i of the subPFs are initialized dependent on xu(j)(t0) and the measurements d^i(t0) for the ith MPC. To initialize the states of xVT,i(j,a)(t0) with a=1,…,Np,j,i of the jth subPF associated with the ith MPC a grid is used. The positions rVT,i(j,a)(t0) of the particles {xVT,i(j,a)(t0)}a=1NP,j,i are distributed on a grid inside a circle around ru(j)(t0) with radius d^i(t0)+Δd such that(29)0≤rVT,ij,at0-rujt0≤d^it0+Δdwith spacing Δd (the number of grid points NP,j,i(tk) can be estimated by Gauss’s circle problem). The additional propagation length is dVT,i(j,a)(t0)=d^i(t0)-rVT,i(j,a)(t0)-ru(j)(t0), where we inherently assume bu(t0)=0 for the initialization. Hence, the total number of particles can be calculated as(30)Nttk=∑j=1Np∑i=0Ntk-1NP,j,itk.For each particle j of the superPF, the receiver state xu(j)(tk-1) is propagated according to the transition model described in Section 3.2 indicated by the function f(xu(j)(tk-1),Ψ˙(tk),nt(tk)) using the heading changes Ψ˙(tk) (cf. Line (5) 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 (7) in Algorithm 1). In case that MPCs are not tracked by KEST anymore, the corresponding subPFs are removed (cf. Line (9) 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 xVT,i(tk) is time-invariant; hence, each subPF assigns the states of the VTs with xVT,i(j,a)(tk)=xVT,i(j,a)(tk-1). Thereafter, the weights for the subPFs and superPF are calculated using (28) and (27).Algorithm 1:Channel-SLAM for time instant tk.
Input:
Multipath estimates: z(tk), σz(tk) and Ψ˙(tk);
States for tk>t0:
x
u
(
j
)
(
t
k
-
1
)
,
x
V
T
,
i
(
j
,
a
)
(
t
k
-
1
)
a
=
1
N
p
,
j
,
i
i
=
0
N
(
t
k
-
1
)
-
1
j
=
1
N
s
Output:
States for ≥t0:
x
u
(
j
)
(
t
k
)
,
x
V
T
,
i
(
j
,
a
)
(
t
k
)
a
=
1
N
p
,
j
,
i
i
=
0
N
(
t
k
)
-
1
j
=
1
N
s
MMSE estimate: x^(tk) for tk>t0
(1) iftk=t0then
(2) Initialization using z(tk) and σz(tk);
(3) else
(4) forj=1:Nsdo
(5) Draw xu(j)(tk)=f(xu(j)(tk-1),Ψ˙(tk),w(tk-1));
(6) ifNew paths detectedthen
(7) Initialize new subPFs;
(8) ifTracking of paths lostthen
(9) Delete corresponding subPFs;
(10) fori=0:N(tk)-1do
(11) fora=1:NP,j,i(tk)do
(12) Assign xVT,i(j,a)(tk)=xVT,i(j,a)(tk-1);
(13) Calculate wij,a(tk)=p(zi(tk)∣xujtk,xVT,ij,atk);
(14) Calculate total subPF weight:
t
j
,
i
=
S
U
M
[
{
w
i
(
j
,
a
)
(
t
k
)
}
a
=
1
N
P
,
j
,
i
(
t
k
)
];
(15) w(j)(tk)=∏i=0N(tk)-1‍tj,i;
(16) Resample using Algorithm 2;
(17) Calculate MMSE x^(tk) according to (31);Algorithm 2:Channel-SLAM resampling algorithm.
Input:
States and weights:
x
u
j
(
t
k
)
,
w
(
j
)
(
t
k
)
,
x
V
T
,
i
(
j
,
a
)
(
t
k
)
,
w
i
(
j
,
a
)
(
t
k
)
a
=
1
N
p
,
j
,
i
(
t
k
)
i
=
0
N
(
t
k
)
-
1
j
=
1
N
s
Output:
Resampled states and weights:
x
~
u
j
(
t
k
)
,
w
~
j
(
t
k
)
,
x
~
V
T
,
i
j
,
a
(
t
k
)
,
w
~
i
(
j
,
a
)
(
t
k
)
a
=
1
N
p
,
j
,
i
(
t
k
)
i
=
0
N
(
t
k
)
-
1
j
=
1
N
s
(1) Initialize the CDF for superPF: cp1=w1(tk);
(2) forj=2:Nsdo
(3) Construct CDF for superPF:
c
p
j
=
c
p
j
-
1
+
w
j
(
t
k
);
(4) f=1;
(5) Draw starting point: up1~U0,Ns-1;
(6) forj=1:Nsdo
(7) upj=up1+Ns-1(j-1);
(8) whileupj>cpfdo
(9) f=f+1;
(10) Assign: x~uj(tk),w~j(tk)=xuftk,1/Ns;
(11) fori=0:N(tk)-1do
(12) Initialize the CDF for (f,i)-th subPF:
c
b
1
=
w
i
(
f
,
1
)
(
t
k
);
(13) fora=2:NP,i,f(tk)do
(14) Construct CDF for subPF:
c
b
a
=
c
b
a
-
1
+
w
i
(
f
,
a
)
(
t
k
);
(15) g=1,ar=1;
(16) fora=1:NP,f,i(tk)do
(17) uba=ub1+1NP,i,f(tk)(a-1);
(18) whileuba>cbgdo
(19) g=g+1;
(20) ifB(x~VT,ij,artk)≤Nmthen
(21) Assign: x~VT,ij,ar(tk)=xVT,if,g(tk);
(22) ar=ar+1;
(23) Bx~VT,ij,artk=Bx~VT,ij,artk+1;
(24) Update number of particles: NP,j,i(tk)=ar;
(25) fora=1:NP,j,i(tk)do
(26) Draw x~VT,ij,atk from the Gaussian-Kernel;
(27) Assign: w~ij,a(tk)=1/NP,j,i(tk);Afterwards, the subPFs and superPF are resampled. The basic idea of the resampling method is to eliminate particles with low weights and reproduce particles with high weights. Algorithm2 shows a pseudocode of the resampling algorithm of Channel-SLAM which is based on the systematic resampling algorithm [49]. Similarly to Algorithm 1, the Channel-SLAM resampling algorithm consists of a resampling algorithm for the superPF which includes resampling algorithms for the subPFs. First, the estimated sampled cumulative distribution function (CDF) of the superPF is constructed, presented by a vector cp with length Np and element cpj with j=1,…,Np. According to the estimated sampled CDF of the superPF, the subPFs are eliminated or resampled. The particles of the superPF with index f are assigned to the resampled particle index j; see Algorithm 2, Line (10), for the assignment of the receiver state. Afterwards, the (f,i)th subPF is resampled with i=0,…,N(tk)-1 using a systematic resampling algorithm, where cb represents the estimated sampled CDF of the (f,i)th subPF.As mentioned before, whenever a new MPC is tracked, many particles are initialized to cover all possible VT positions in each subPF. For example if theith MPC has a delay of d^i(t0) = 30 m, each ith subPF would use NP,j,i(t0)=2821 particles according to (29) for j=1,…,Np with spacing Δd = 1 m. However during the receiver movement many particles of the subPFs are resampled at the same grid point because the states of the VTs xVT,i(tk) are time-invariant. In order to adapt the number of particles, we limit the number of resampled particles per grid point to Nm; see Algorithm 2, Line (20). Evaluations in Section 5 show that the reduction of the number of particles does not influence the positioning accuracy but, however, leads to a gain on computational performance. Afterwards, the states of the VTs xVT,i(j,a)(tk) are drawn using a Gaussian-Kernel (cf. Line (26) in Algorithm 2). As stated in [31], the Gaussian-Kernel has a low bandwidth which does not influence the grid based reduction method.Usually, we are interested in a point estimate of the state instead of its a posteriori PDF. According to [31], the MMSE point estimate, x^(tk)=x^u(tk)T,x^VT(tk)TT (cf. Line (17) in Algorithm 1), is calculated by(31)x^tk≈∑j=1Nswjtk×xujtk∑a=1NP,j,0tkw0j,atkxVT,0j,atk⋮∑a=1NP,j,Ntk-1tkwNtk-1j,atkxVT,Ntk-1j,atk.


## 4. Posterior Cramér-Rao Lower Bound

The PCRLB can be calculated by the inverse of the posterior information matrixJ(tk) and provides a lower bound of the variance of a Bayesian estimator [50] with(32)Ex^tk-xtkx^tk-xtkT=Mtk≥Jtk-1.The inequality in (32) means that the difference M(tk)-J(tk)-1 is a positive semidefinite matrix. For the performance evaluation of a filter like Channel-SLAM with the system equations of (10) and (11), the posterior information matrix can be calculated recursively according to [51], with(33)Jtk=D22tk-D21tkJtk-1+D11tk-1D12tk,where(34)D11tk=-EΔxtk-1xtk-1ln⁡pxtk∣xtk-1,Ψ˙tk,D21tk=-EΔxtkxtk-1ln⁡pxtk∣xtk-1,Ψ˙tk=D12tkT,D22tk=-EΔxtkxtkln⁡pxtk∣xtk-1,Ψ˙tk-EΔxtkxtkln⁡pztk∣xtk,where ∇a stands for the first-order partial derivatives with respect to a and Δab stands for the second-order partial derivatives with Δab≜∇a∇bT. In order to calculate the PCRLB, we use the transition model of Section 3.2.1. In case of non-Gaussian noise and nonlinearity as in Section 3.2.2, the expectation estimator in (34) has to be approximated by Monte Carlo simulations. To calculate the PCRLB, we combine the time-invariant transition model for the VTs xVT(tk) as introduced in (14) and the transition model of the receiver (35), with(35)xtk=Autδ,Ψ˙tk00I︸Atδ,Ψ˙tkxtk-1+nttk,where nt(tk)~N0,Q(tδ) is the transition noise with covariance matrix(36)Qtδ=Qutδ000,where Qu(tδ) is defined in (21).Under the condition of a known control inputΨ˙(tk), see [51], and using the linear transition model of Section 3.2.1 and Gaussian distributed transition noise of (36), we obtain for (34)(37)D11tk=Atδ,Ψ˙tkTQtδ-1Atδ,Ψ˙tk,D12tk=-Atδ,Ψ˙tkTQtδ-1D22tk=Qtδ-1+EFtk,where the matrix F(tk) is the snapshot based Fisher information matrix. Substituting (37) into (33), we obtain(38)Jtk=EFtk+Qtδ+Atδ,Ψ˙tkJtk-1-1Atδ,Ψ˙tkT-1,using the matrix inversion lemma because of the singularity of Q(tδ).The snapshot based Fisher information matrixF(tk) in (38) can be obtained by(39)Ftkk,w=∂μxtkH∂xtkkR-1tk∂μxtk∂xtkw,with the covariance matrix R(tk), and(40)μxtk=d0tk,…,dNtk-1tk,with the propagation lengths di(tk) of (16) for the ith MPC.


## 5. Evaluations Based on Measurements

This section evaluates the derived algorithm based on channel measurements on an airfield in front of a hangar with a single static physical transmitter and a moving pedestrian as shown in Figure6. Figure 6 provides the scenario by a top view with the physical transmitter position in red, the track in blue, the starting position in green, and the end position in magenta. The measurements are conducted using the MEDAV RUSK-DLR broadband channel sounder in single-input single-output (SISO) mode with the measurement setup as indicated in Figure 7. The transmitter and receiver are connected to the same rubidium clock to prevent time drifts during the measurements. The static physical transmitter emits a 10 mW multitone signal (see [52]) with N=1281 subcarriers having equal gains at a center frequency of 1.51 GHz and a bandwidth of B=100MHz. On the receiver side, the CIR snapshots are repeatedly measured every Tg=1.024ms. As shown in Figure 7, the receiving antenna is mounted on a pole attached on the backpack of the pedestrian. Additionally, the pedestrian is equipped with a hand-held equipment including a Xsense IMU (MTI-G-700) and a laptop which stores the IMU measurements. In order to obtain the ground truth of the pedestrians movement, a prism is mounted next to the antenna at the pole above the pedestrian. The prism is tracked by a tachymeter (TPS1200 from Leica Geosystems AG) which sends the measured coordinates to the laptop that records the coordinates simultaneously with the IMU measurements. The tachymeter gives a nominal accuracy in the subcentimeter domain. To synchronize all devices, the laptop is additionally connected by cable to the channel sounder. Thus, we are able to obtain the ground truth of the pedestrian for each captured CIR snapshot. Although the synchronization between the IMU and the channel sounder might be in the ms scale only, the influence on the position estimation is negligible because of the low pedestrian speed of around 0.7 m/s.Figure 6
Measurement scenario with a fixed transmitter and a moving receiver (pedestrian). The pedestrian moves on the blue track for 155 s, in total 111 m. The starting position and end position are indicated by the green and magenta circles. The metalized doors of the hangar and the chain-link fences act as reflecting surfaces for the transmitted wireless signal. Hence, we can calculate the corresponding VT positions by mirroring the physical transmitter position on the reflecting surfaces.Figure 7
Measurement setup: transmitter and receiver use the same rubidium clock for synchronization. The receiving antenna is mounted on a stick next to a prism for measuring the ground truth of the moving pedestrian. Additionally, the pedestrian is holding a hand-held device which consists of an IMU and a laptop.The pedestrian is moving on the indicated blue track of Figure6 for 155 s or 111 m in front of a hangar with metalized doors. During the whole pedestrian movement, the LoS path between transmitter and receiver is present. Figure 8 shows the recorded CIRs versus the pedestrian moving time in seconds, where the time delays of the CIRs are multiplied by the speed of light, thus, in meters.Figure 8
Recorded CIRs versus the pedestrian moving time in seconds.In order to exploit the multipath propagation for positioning, we have to estimate and track the MPCs over time. Hence, the accuracy of Channel-SLAM relies directly on the accuracy of the CIR estimations of KEST. Channel-SLAM considers an underdetermined system; therefore, long visible paths are preferable. Thus for the evaluations, we extract from the KEST estimates only the long visible paths as visualized in Figure9 (Channel-SLAM could use all detected MPCs; however, this would increase the computational complexity). Figure 9 shows the estimation results of KEST for the CIR versus the pedestrian moving time in seconds. The black circled line in Figure 9 indicates the geometrical line-of-sight (GLoS) path delay, which matches perfectly with the KEST estimates for the first path. Additionally, other MPCs can be tracked by KEST for a long time. The metalized doors of the hangar act as a reflecting surface for the transmitted wireless signal. We can obtain the position of VT1 by mirroring the physical transmitter on the reflecting surface as mentioned in Section 2. Additionally, the chain-link fences indicated by Fence 1, Fence 2, and Fence 3 act as reflecting surfaces. We obtain VT2, VT3, and VT4 by mirroring the physical transmitter position at Fence 1, Fence 2, and Fence 3. The positions of the hangar, Fence 1, Fence 2, and Fence 3 are measured using the tachymeter; thus, we are able to calculate the positions of VT1, VT2, VT3, and VT4. Please note that VT4 is not shown in Figure 6. Based on the calculated VT positions, we are able to calculate the hypothetical propagation distances between these VTs and the moving pedestrian. We can see that they match the KEST estimates as indicated by the black lines in Figure 9. The measurement scenario considers only one time signal reflections. For examples on VTs with multiple number of reflections, diffractions, or scattering, please see [31].Figure 9
Estimation results of KEST for the CIR versus the pedestrian moving time in seconds where only long tracked paths are visualized. Additionally it shows the calculated propagation paths of the GLoS path and paths ofVT1, VT2, VT3, and VT4.The evaluations are performed usingNs=2000 particles in the superPF, whereas the number of particles for the subPFs for each MPC is different depending on the estimated delay of each MPC. Channel-SLAM obtains the measurements z(tk) from KEST and the heading change Ψ˙(tk) from the IMU every tδ=0.1s. For the initialization of Channel-SLAM, we use prior information p(xu(t0)). The prior information includes the starting position and moving direction, whereas the speed is initialized using a uniform distribution between 0 m/s and 1 m/s. Please note that an unknown starting position and direction or larger initial uncertainties may result in a biased and rotated coordinate system in the estimation. We empirically set Δd = 1 m. Since Channel-SLAM has no knowledge of the physical transmitter position, Channel-SLAM estimates the position of the physical transmitter as a VT. During the pedestrian movement, the number of tracked MPCs changes which results in removing and initialization of subPFs during the movement. Additionally, because Channel-SLAM does not consider retracking of previous MPCs, for example, the MPCs of VT1 and VT3 which are tracked multiple times, they are initialized without any prior information.To see the positioning performance of Channel-SLAM, we compare the following algorithms and bounds:(i)
Dynamic-Channel-SLAM: it is a Channel-SLAM implementation using the dynamical adaptation of the number of particles as introduced in Section3.4 where we limit the number of particles per bin to Nm=30 and the grid size to Δd=1m.(ii)
RBPF-Channel-SLAM: it is similar to Dynamic-Channel-SLAM, however, without using the dynamical adaptation of the number of particles.(iii)
VT-Knowledge-Algo: it is a positioning algorithm with perfect knowledge of all VT positionsrVT,i(tk) and additional propagation lengths dVT,i. Because the measurement scenario considers only one time reflections, VT-Knowledge-Algo reflects algorithms of [24, 53] which consider reflected signals as signals emitted from VTs, where the VT positions are precalculated based on the knowledge of the reflecting surface and physical transmitter positions. VT-Knowledge-Algo can be seen as a lower bound for Channel-SLAM. Similarly to Channel-SLAM, VT-Knowledge-Algo uses the delays of the estimated MPCs provided by KEST as input, assumes the knowledge of starting position and direction, and is implemented using a PF with Ns=2000 particles.(iv)
PCRLB: it is as the PCRLB derived in Section4 using σqu2=5·10-4m2/s3, the standard deviation σz(tk) for each MPC from KEST, and the prior like the Channel-SLAM algorithms.The above-mentioned algorithms are evaluated using the Gaussian-Transition-Model of Section3.2.1 indicated by index 1 and the Rician-Transition-Model of Section 3.2.2 indicated by index 2, for example, RBPF-Channel-SLAM-1 and RBPF-Channel-SLAM-2 for RBPF-Channel-SLAM. For the Gaussian-Transition-Model we set the continuous-time process noise intensity to σqu2=5·10-4m2/s3 and for the Rician-Transition-Model we set the standard deviation of the heading noise of (25) to σΦ(tk)=0.1° and the scale parameter of the velocity of (24) to σm(tk)=0.025m/s.Figure10 shows the root mean square error (RMSE) RMSEu(tk)=E{rutk-r^utk2} of the estimated pedestrian position r^u(tk) versus the pedestrian moving time for Dynamic-Channel-SLAM-1 in magenta, RBPF-Channel-SLAM-1 in cyan, and VT-Knowledge-Algo-1 in yellow and the black line indicates the PCRLB. Because the PF includes randomness, the position estimates differ for each evaluation unless the number of particles is infinite even if the same measurement data are used. Therefore, we perform 200 independent evaluations based on the same measurement data. For the evaluations we add an artificial clock bias to the measurements to verify the clock bias estimation capabilities. Because of the initialization of the receiver position using prior knowledge, all algorithms perform similarly at the beginning of the track where the position error is rather low. Afterwards, RMSEu(tk) is varying between 0.6 m and 4 m for Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1. VT-Knowledge-Algo-1 can be interpreted as a lower bound and estimates the receiver position with the lowest RMSE. Between 70 s and 90 s of the receiver movement, VT-Knowledge-Algo-1 has a slightly higher RMSE which might be due to the nonperfect reflecting surfaces, KEST estimation errors, or small inaccuracies in the calculations of the VT positions. Furthermore, we see that we obtain a similar RMSE for Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1. However, if we have a look on the number of used particles, as shown in Figure 11, we see a major computational performance gain of the Dynamic-Channel-SLAM-1 compared to RBPF-Channel-SLAM-1. Figure 11 shows the total number of particles Nt(tk) calculated according to (30) for Np=2000 versus the pedestrian moving time. At the beginning, both Channel-SLAM algorithms are initialized with the same number of particles. As soon as the pedestrian is moving, the estimations of the VT positions are converging resulting in a reduction of the number of particles for Dynamic-Channel-SLAM-1. When new MPCs are tracked (see Figures 9 and 11), the total number of particles Nt(tk) increases and reduces afterwards during runtime. Especially at the end of the track, Dynamic-Channel-SLAM-1 uses 40 times less particles than RBPF-Channel-SLAM-1.Figure 10
R
M
S
E
u
(
t
k
) versus the pedestrian moving time in seconds for Dynamic-Channel-SLAM-1, RBPF-Channel-SLAM-1, VT-Knowledge-Algo-1, and the PCRLB.Figure 11
Total number of particlesNt(tk) versus the pedestrian moving time in seconds for Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1.As mentioned before, the black line in Figure10 indicates the PCRLB. The PCRLB shows the theoretical performance bound which has on average a 2-3 times lower RMSE than Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1. However, the curve shapes of the PCRLB, Dynamic-Channel-SLAM-1, and RBPF-Channel-SLAM-1 are similar. By increasing the number of particles, the RMSE of Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1 might be decreased. Additionally, the PCRLB shows a theoretical limit which is not affected by estimation inaccuracies of KEST caused by nonperfect reflecting surfaces or DMCs.Figure12 shows the RMSE for Dynamic-Channel-SLAM-1 in magenta, RBPF-Channel-SLAM-1 in cyan, VT-Knowledge-Algo-1 in yellow, Dynamic-Channel-SLAM-2 in blue, RBPF-Channel-SLAM-2 in red, and VT-Knowledge-Algo-2 in green. Similarly to Figure 12, we see that we obtain a similar RMSE for Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1. Additionally, we see a significant performance gain in the position accuracy by comparing the different transition models. Similarly to Figure 12, VT-Knowledge-Algo-1 has between 80 s and 90 s a slightly higher RMSE because of the same reasons stated before. Furthermore, we see as well that we obtain a similar RMSE for Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1.Figure 12
R
M
S
E
u
(
t
k
) versus the pedestrian moving time in seconds for Dynamic-Channel-SLAM-1, RBPF-Channel-SLAM-1, VT-Knowledge-Algo-1, Dynamic-Channel-SLAM-2, RBPF-Channel-SLAM-2, and VT-Knowledge-Algo-2.Figure13 shows the enlarged measurement scenario with estimated PDFs for the physical transmitter, VT1, VT2, and the pedestrian position. Whereas the PDFs of the physical transmitter, VT2, and pedestrian position are the estimation results at the end of the track, the PDFs of VT1 are the estimation results when the tracking of the MPC is lost, that is, after 75 s. We see that especially the physical transmitter and VT2 position can be estimated with a low uncertainty. Additionally, Figure 13 shows two examples of the MMSE point estimates of the receiver position for Dynamic-Channel-SLAM-2 in red, VT-Knowledge-Algo-2 in green, and ground truth of the track in blue. We see that we obtain similar position estimation results for both algorithms.Figure 13
Enlarged measurement scenario of Figure6 with the ground truth of the track in blue. The figure shows two examples of the MMSE point estimates of the receiver position for Dynamic-Channel-SLAM-2 in red and VT-Knowledge-Algo-2 in green. The starting position and end position are indicated by the green and magenta circles. Additionally, it shows the estimated PDFs for the physical transmitter, VT1, VT2, and pedestrian position.


## 6. Conclusions

In this paper, we extended the work on multipath assisted positioning, called Channel-SLAM. The new positioning method takes advantage of the multipath components (MPCs) instead of mitigating them. In the proposed approach, multipath signals are treated as signals from virtual transmitters (VTs), where the locations of these VTs are unknown. To improve the accuracy of Channel-SLAM, an inertial measurement unit (IMU) is used to obtain heading information of the moving receiver. Furthermore, we investigate a novel particle filter (PF) implementation which adapts the number of particles during runtime. Measurement results show that accurate position estimation is possible without the knowledge of the physical transmitter position. Hence, the new algorithm does not rely on prior information such as the room layout or information collected in a database for fingerprinting except for the initial position and direction. We compare the position accuracy of Channel-SLAM to that of a derived posterior Cramér-Rao lower bound (PCRLB). Additionally, we obtain similar position estimation results with Channel-SLAM and an algorithm with perfect knowledge of all VT positions.



---

*Source: 9170746-2017-10-02.xml*