# Rao-Blackwellized Gaussian Sum Particle Filtering for Multipath Assisted Positioning

Authors: Markus Ulmschneider; Christian Gentner; Thomas Jost; Armin Dammann
Journal: Journal of Electrical and Computer Engineering
 (2018)

---

## Abstract

In multipath assisted positioning, multipath components arriving at a receiver are regarded as being transmitted by a virtual transmitter in a line-of-sight condition. As the locations and clock offsets of the virtual and physical transmitters are in general unknown, simultaneous localization and mapping (SLAM) schemes can be applied to simultaneously localize a user and estimate the states of physical and virtual transmitters as landmarks. Hence, multipath assisted positioning enables localizing a user with only one physical transmitter depending on the scenario. In this paper, we present and derive a novel filtering approach for our multipath assisted positioning algorithm called Channel-SLAM. Making use of Rao-Blackwellization, the location of a user is tracked by a particle filter, and each landmark is represented by a sum of Gaussian probability density functions, whose parameters are estimated by unscented Kalman filters. Since data association, that is, finding correspondences among landmarks, is essential for robust long-term SLAM, we also derive a data association scheme. We evaluate our filtering approach for multipath assisted positioning by simulations in an urban scenario and by outdoor measurements.


---

## Body

## 1. Introduction

The amount of available and potential services requiring precise localization of a user has steadily increased over the recent years. Global navigation satellite systems (GNSSs) can often satisfy the demands for localization in scenarios where the receiver has a clear view of the sky. However, if the view of the sky is obstructed, such as indoors, in urban canyons, or in tunnels, the positioning performance of GNSSs may be drastically decreased, or no positioning solution may be obtained at all [1]. Reasons for this include a low received signal power due to signal blocking or shadowing and multipath propagation.In contrast to GNSS signals, many kinds of terrestrial signals are likely to have a good coverage in GNSS denied places. In particular, cellular radio frequency (RF) signals are designed to be reliably available at least in populated areas, and they may be used as signals of opportunity (SoOs) for positioning. However, also terrestrial signals experience multipath propagation. Multipath propagation biases range estimates if standard correlator based methods are used. Various approaches to handle the multipath problem have been addressed in the literature, for example, in [2]. Advanced methods such as maximum likelihood (ML) mitigation algorithms try to estimate the channel impulse response (CIR) and to mitigate the influence of multipath components (MPCs) on the line-of-sight (LoS) path [3].The idea of multipath assisted positioning is contrary, though. Instead of regarding multipath propagation as ill, the spatial information of MPCs on the receiver position is exploited. In [4], the information of MPCs is used in a fingerprinting scheme. Going one step further, each MPC can be regarded as being transmitted by a virtual transmitter in a pure LoS condition, and the virtual transmitters can be used to locate the user. Such an approach is called multipath assisted positioning.The authors of [5, 6] derived some theoretical bounds for multipath assisted positioning. Multipath assisted positioning schemes have, for example, been applied in radar applications [7], using ultrawideband (UWB) [8, 9] or 5G [10] systems and in cooperative systems [11].If the locations of physical transmitters and reflecting and scattering objects are known, the locations of virtual transmitters can be calculated based on geometrical considerations. The authors of [12] assume the room layout to be known and focus on the association among virtual transmitters and reflecting walls. In a general setting, however, the scenario is unknown to the user.The authors of [13, 14] have presented a multipath assisted positioning scheme named Channel-SLAM that does not rely on prior information on the scenario. Instead, the locations of the physical and virtual transmitters are estimated simultaneously with the user position in a simultaneous localization and mapping (SLAM) [15, 16] approach. In general, SLAM describes the simultaneous estimation of a user position and the locations of landmarks. In Channel-SLAM, the landmarks are the physical and virtual transmitters. Previous extensions to Channel-SLAM include mapping of the user positions [17], the consideration of vehicular applications [18], and data association methods [19, 20], for example.Nonlinearities in the prediction and update equations of the Bayesian recursive estimation framework prohibit the use of optimal algorithms such as the Kalman filter, since the integrals involved in the estimation process cannot be solved in closed form or become intractable. A popular alternative is the extended Kalman filter (EKF) [21], which linearizes the nonlinear terms using a first-order Taylor series expansion. However, such a linearization can introduce large errors in the estimation process [22]. The unscented Kalman filter (UKF) [23, 24] uses a nonlinear transformation to deal with nonlinearities and outperforms the EKF in a wide range of applications [22, 25].UKF methods have found their way into localization problems, for example, in [27, 28]. The authors of [29] propose Gaussian sum cubature filters. In [30, 31], the authors consider a Rao-Blackwellization scheme for SLAM with a particle filter for the user state and UKFs for the landmark states, where the measurement model is based on linearization, though.The current Channel-SLAM algorithm uses a Rao-Blackwellized particle filter to estimate the user state and the location of transmitters simultaneously. Hence, both the user state probability density function (PDF) and the transmitter state PDFs are represented by a large set of particles, tending to result in a high memory occupation. This paper is an extension of [32], where we proposed a novel estimation approach for Channel-SLAM scheme based on Rao-Blackwellization and performed first simulations. We refer to this new estimation method as Rao-Blackwellized Gaussian sum particle filter (RBGSPF). In the RBGSPF, the user position is tracked by a sequential importance resampling (SIR) particle filter, while the physical and virtual transmitter state PDFs are represented by Gaussian mixture models estimated by UKFs. This parametrized representation of the transmitter states is a key enabler for exchanging maps of transmitters among users, since the amount of data that has to be communicated among users can be decreased drastically compared to the nonparametric representation with particles. Such an exchange of maps may be performed directly among users or via a central entity, for example, in form of local dynamic maps (LDMs) in an intelligent transportation system (ITS) context. In this paper, we provide a full and detailed derivation of our novel algorithm. In particular, we derive the calculation of the particle weights in the user particle filter given the representation of the transmitters in the UKF framework. Since data association is an essential feature for the accuracy in long-term SLAM, we also derive a data association method based on [33]. We evaluate our algorithm by both simulations in an urban scenario and outdoor measurements.The remainder of this paper is structured as follows. Section2 describes the fundamental idea behind multipath assisted positioning and Channel-SLAM. In Section 3, we briefly summarize some concepts of nonlinear Kalman filtering. The derivation of the RBGSPF is presented in Section 4, and a solution to data association is presented in Section 5. After the experimental results in Section 6, Section 7 concludes the paper.Throughout the paper, we use the following notation:(i)
As indices,i stands for a user particle, j denotes a transmitter or a signal component, l is a component in a Gaussian mixture model, and m stands for a sigma point.(ii)
(·)T denotes the transpose of a matrix or vector.(iii)
1n denotes the identity matrix of dimension n×n.(iv)
0n and 0m×n denote the zero matrices of dimensions n×n and m×n, respectively.(v)
N(x;μ,C) denotes the PDF of a normal distribution in x with mean μ and covariance C.(vi)
c0 denotes the speed of light.(vii)
· denotes the Euclidean norm of a vector.


## 2. Multipath Assisted Positioning

### 2.1. Virtual Transmitters

The idea of virtual transmitters is illustrated in Figure1. The physical transmitter Tx transmits an RF signal. A mobile user equipped with an RF receiver receives the transmitted signal via three different propagation paths.Figure 1
Signals from the physical transmitter Tx are received at the two user positions via different propagation paths. Each MPC arriving at the receiver is regarded as being transmitted by a virtual transmitter in a pure LoS condition. The propagation paths correspond to a reflection at the wall (vTx1), a scattering at a point scatterer (vTx2), and a scattering followed by a reflection at the wall (vTx3).In the first case, the signal is reflected at the reflecting surface. The user treats the corresponding impinging MPC as being transmitted by the virtual transmitter vTx1 in a pure LoS condition. The location of vTx1 is the location of the physical transmitter Tx mirrored at the reflecting surface. When the user moves along the trajectory, the reflection point at the wall moves as well. However, the location of vTx1 is static. The two transmitters Tx and vTx1 are inherently perfectly synchronized.In the second case, the signal from the physical transmitter is scattered by a point scatterer and then received by the user. We define the effect of scattering such that the energy of an electromagnetic wave impinging against the scatterer is distributed uniformly in all directions [34]. The user regards the scattered signal as a LoS signal from the virtual transmitter vTx2, which is located at the scatterer location. If the signal is scattered, the physical and the virtual transmitter are not time synchronized: the virtual transmitter has an additional delay offset to the physical transmitter corresponding to the propagation time of the signal traveling from the physical to the virtual transmitter.In the third case, the signal is first scattered at the scatterer and then reflected at the surface. The user treats this signal as being sent from the virtual transmitter vTx3. The location of vTx3 is the location of vTx2, that is, the scatterer location, mirrored at the reflecting surface. Accordingly, the concept of single reflections and scatterings can be generalized in a straightforward manner to the case of multiple reflections and scatterings by applying the first two cases iteratively. In case the signal undergoes only reflections, the physical and the virtual transmitters are inherently time synchronized. If the signal is scattered at least once, the delay offset corresponds to the actual propagation time of the signal from the physical transmitter to the last scatterer the signal interacts with. Therefore, in Figure1, the virtual transmitters vTx2 and vTx3 have the same delay offset towards the physical transmitter. Note that a delay offset can be interpreted as a clock offset.Throughout the paper, we consider the physical transmitter and the environment to be static. Hence, the virtual transmitters are static as well.


### 2.2. Recursive Bayesian Estimation

Recursive Bayesian estimation [35] is a method to recursively estimate the evolution of a state vector x, where the state evolution is modeled as(1)xk=fkxk-1,vk-1.The index k denotes the time instant, the function fk(·) is assumed to be known, and vk-1 denotes a sample of the process noise with covariance matrix Q. The state is related to the measurement zk by(2)zk=hkxk,nk,where hk(·) is again a known function and nk is a sample of the measurement noise with covariance matrix R. Recursive Bayesian estimation works in two steps, the prediction and the update step. The corresponding PDFs can be calculated recursively by(3)pxk∣z1:k-1=∫pxk∣xk-1pxk-1∣z1:k-1dxk-1for the prediction step and by(4)pxk∣z1:k=1ckpzk∣xkpxk∣z1:k-1,for the update step, where ck is a constant and z1:k denotes the measurements from time instant 1 to k. The state transition prior p(xk∣xk-1) and the measurement likelihood p(zk∣xk) are obtained from the movement model in (1) and the measurement model in (2), respectively.


### 2.3. Channel-SLAM

In the following, we will revise the Channel-SLAM algorithm from [13, 17]. Figure 2 gives an overview of the two stages of Channel-SLAM. In the first stage, the parameters of the signal components received by the user via different propagation paths are estimated. The resulting estimates are used as measurement input in the second stage, where the states of the physical and virtual transmitters and the user position are estimated simultaneously in a SLAM scheme. Further sensors, such as an inertial measurement unit (IMU), may be included in the second stage. The locations of both the physical and virtual transmitters are assumed to be unknown. Thus, Channel-SLAM does not differentiate between physical and virtual transmitters, and the term transmitter comprises both physical and virtual transmitters in the following. Each signal component arriving at the receiver corresponds to one transmitter.Figure 2
Based on the received signal, the parameters of the propagation paths are estimated in the first step by the KEST algorithm. In the second step, the estimates serve as measurements for estimating the positions of the user and the physical and virtual transmitters. In addition, user heading rates of change measurements from an IMU are incorporated in the second step.The RF propagation channel between the physical transmitter and the user equipped with a receiver is assumed to be a linear and time-variant multipath channel.The received signal is modeled as a superposition of signal components of the transmit signals(t), where the jth signal component is defined by a complex amplitude aj(tk) and a delay dj(tk) at time tk. The signal received by the user at time instant tk is(5)yτ,tk=∑jajtksτ-djtk+nτ,where n(τ) is a sample from a colored noise sequence incorporating both dense multipath components (DMCs) and additive Gaussian noise. The channel is assumed to be constant during the short time interval when the received signal is sampled at time instant k.The physical transmitter continuously broadcasts the signals(t) that is known to the user. At the user side, the parameters of the signal components arriving at the receiver are estimated. Such parameters can in general be the complex amplitude, time of arrival (ToA), angle of arrival (AoA), or Doppler shift, depending on the available hardware and the scenario. For the signal parameter estimation, we use the KEST algorithm [36]. The KEST estimator works in two stages. In an inner stage, a ML parameter estimator such as Space-Alternating Generalized Expectation-Maximization (SAGE) [37] estimates the parameters of the signal components jointly on a snapshot basis. An outer stage tracks these estimated parameters over time with a Kalman filter and keeps track of the number of signal components. The KEST estimator is in general able to handle the DMCs in the noise term in (5). However, DMC handling is not implemented in our evaluations, leading to a model mismatch in KEST and hence to a higher variance in the parameter estimation. However, we do not expect many DMCs in our evaluation scenarios. In an indoor scenario, for example, DMCs need to be considered [38].In the literature, there are alternatives to the KEST estimator. For example, the authors of [39] track signal component parameters based on an EKF, though the authors of [40] showed that the KEST estimator is more robust in resolving signal components that are close to each other in the state space. In [41], an EKF is used as well for parameter estimation, while the position estimation is based on the time difference of arrival (TDoA) of virtual transmitters. The authors of [42] consider the linearization of the observation model in the EKF a major drawback that might lead to a tracking loss.In the second stage of Channel-SLAM, we use only the delays, that is, ToAs, and AoAs, estimated by KEST as measurement inputs. Hence, after sampling the received signal, the KEST estimates at time instantk are comprised in the vector(6)zk=dkTθkTT,where(7)dk=d1,k⋯dNTX,kTare the ToA estimates for the NTX signal components, or transmitters, and(8)θk=θ1,k⋯θNTX,kTare the corresponding AoA estimates. Note that the number of signal components and thus transmitters may change over time. Nevertheless, for notational convenience, we omit the time instant index k in NTX.In the second stage of Channel-SLAM, the user statexu,k is estimated simultaneously with the state of the transmitters xTX,k. The entire state vector is hence(9)xk=xu,kTxTX,kTT=xu,kTxTX,k1T⋯xTX,kNTXTT,where xTX,kj is the state of the jth transmitter. As we consider a two-dimensional scenario, the user state at time instant k is defined by(10)xu,k=xkykvx,kvy,kT=pu,kTvu,kTT,where the user position is defined by pu,k=[xkyk]T and the user velocity by vu,k=[vx,kvy,k]T. Each transmitter is defined by its location pTX,k=[xTX,kyTX,k]T and a clock offset τ0,k at time instant k. The state vector of the jth transmitter is hence defined by(11)xTX,kj=xTX,kjyTX,kjτ0,kjT=pTX,kjTτ0,kjT.Our goal is to find the minimum mean square error (MMSE) estimator forxk, which is defined as (12)x^k=∫xkpxk∣z1:kdxk, where z1:k denotes all measurements up to time instant k.We use a recursive Bayesian estimation scheme as in Section2.2 to estimate the posterior PDF p(xk∣z1:k). This posterior can be factorized as(13)pxk∣z1:k=pxTX,k,xu,k∣z1:k=pxTX,k∣xu,k,z1:kpxu,k∣z1:k.The signal components arriving at the receiver are assumed to be independent of each other; that is, we assume they interact with distinct objects. Assuming independence among the measurements for distinct transmitters (on the one hand, the parameters of the signal components are estimated jointly by the KEST algorithm, and hence these estimates might be correlated between signal components and between the parameters; on the other hand, the correlation is likely to have effect only on a short term basis as KEST estimates are unbiased when observed over a longer time), the first factor in (13) can be factorized further as(14)pxTX,0:k∣xu,0:k,z1:k=∏j=1NTXpxTX,0:kj∣xu,0:k,z1:kj.With the above factorization, the transmitter states are estimated independently from each other.As we consider a static scenario, the virtual transmitters are static as well, and the transition prior for thejth transmitter is calculated as(15)pxTX,kj∣xTX,k-1j=δxTX,kj-xTX,k-1j,where δ(·) denotes the Dirac distribution.For the prediction of the user, additional sensors such as an IMU carried by the user may be integrated into the movement model. Within this paper, we assume only heading change rate measurements from a gyroscope to be available, though, and no knowledge on the user speed. With the gyroscope heading change rateΔβ,k, we predict the movement of the user by(16)xu,k=12Tk1202RΔβ,kxu,k-1=Fu,kxu,k-1,where Tk denotes the time between instants k-1 and k. The two-dimensional rotation matrix R(Δβ,k) is defined as(17)RΔβ,k=cos⁡Δβ,k+wk-sin⁡Δβ,k+wksin⁡Δβ,k+wkcos⁡Δβ,k+wk,where wk is the heading noise which is distributed following a von Mises distribution. Hence, the function fk in (1) can be expressed in our case as(18)fkxk-1,vk-1=Fu,k04×3NTX03NTX×413NTXxk-1+vk-1,where the process noise covariance matrix Q is diagonal.As depicted in Figure3, an AoA measurement for a transmitter j describes the angle θj,k between the user heading direction vu,k and the incoming signal from the transmitter. The measurement noise for the ToA and AoA measurements is assumed to be zero-mean Gaussian distributed with variances σd,j2 and σθ,j2, respectively, for the jth transmitter. Also, we assume no cross-correlation between the single ToA and AoA measurements. The likelihood for the measurement vector zk conditioned on the state vector xk is therefore the product(19)pzk∣xk=∏j=1NTXNdj,k;d^j,k,σd,j2Nθj,k;θ^j,k,σθ,j2,where the predicted ToA between the user and the jth transmitter is(20)d^j,k=1c0pu,k-pTX,kj+τ0,kj,and the predicted AoA is calculated as(21)θ^j,k=atan2⁡yk-yTX,kj,xk-xTX,kj-atan2⁡vy,k,vx,k.The function atan2(y,x) calculates the four-quadrant inverse tangent function. It returns the counterclockwise angle between the positive x-axis and the point given by the coordinates (x,y).Figure 3
The user moves in the directionvu,k at time instant k. The heading change rate from the IMU is Δβ,k. The ToA and AoA measurements for the signal from the jth transmitter are dj,k and θj,k, respectively, where θj,k describes the angle between the heading direction vu,k of the user and the arriving signal.


## 2.1. Virtual Transmitters

The idea of virtual transmitters is illustrated in Figure1. The physical transmitter Tx transmits an RF signal. A mobile user equipped with an RF receiver receives the transmitted signal via three different propagation paths.Figure 1
Signals from the physical transmitter Tx are received at the two user positions via different propagation paths. Each MPC arriving at the receiver is regarded as being transmitted by a virtual transmitter in a pure LoS condition. The propagation paths correspond to a reflection at the wall (vTx1), a scattering at a point scatterer (vTx2), and a scattering followed by a reflection at the wall (vTx3).In the first case, the signal is reflected at the reflecting surface. The user treats the corresponding impinging MPC as being transmitted by the virtual transmitter vTx1 in a pure LoS condition. The location of vTx1 is the location of the physical transmitter Tx mirrored at the reflecting surface. When the user moves along the trajectory, the reflection point at the wall moves as well. However, the location of vTx1 is static. The two transmitters Tx and vTx1 are inherently perfectly synchronized.In the second case, the signal from the physical transmitter is scattered by a point scatterer and then received by the user. We define the effect of scattering such that the energy of an electromagnetic wave impinging against the scatterer is distributed uniformly in all directions [34]. The user regards the scattered signal as a LoS signal from the virtual transmitter vTx2, which is located at the scatterer location. If the signal is scattered, the physical and the virtual transmitter are not time synchronized: the virtual transmitter has an additional delay offset to the physical transmitter corresponding to the propagation time of the signal traveling from the physical to the virtual transmitter.In the third case, the signal is first scattered at the scatterer and then reflected at the surface. The user treats this signal as being sent from the virtual transmitter vTx3. The location of vTx3 is the location of vTx2, that is, the scatterer location, mirrored at the reflecting surface. Accordingly, the concept of single reflections and scatterings can be generalized in a straightforward manner to the case of multiple reflections and scatterings by applying the first two cases iteratively. In case the signal undergoes only reflections, the physical and the virtual transmitters are inherently time synchronized. If the signal is scattered at least once, the delay offset corresponds to the actual propagation time of the signal from the physical transmitter to the last scatterer the signal interacts with. Therefore, in Figure1, the virtual transmitters vTx2 and vTx3 have the same delay offset towards the physical transmitter. Note that a delay offset can be interpreted as a clock offset.Throughout the paper, we consider the physical transmitter and the environment to be static. Hence, the virtual transmitters are static as well.


## 2.2. Recursive Bayesian Estimation

Recursive Bayesian estimation [35] is a method to recursively estimate the evolution of a state vector x, where the state evolution is modeled as(1)xk=fkxk-1,vk-1.The index k denotes the time instant, the function fk(·) is assumed to be known, and vk-1 denotes a sample of the process noise with covariance matrix Q. The state is related to the measurement zk by(2)zk=hkxk,nk,where hk(·) is again a known function and nk is a sample of the measurement noise with covariance matrix R. Recursive Bayesian estimation works in two steps, the prediction and the update step. The corresponding PDFs can be calculated recursively by(3)pxk∣z1:k-1=∫pxk∣xk-1pxk-1∣z1:k-1dxk-1for the prediction step and by(4)pxk∣z1:k=1ckpzk∣xkpxk∣z1:k-1,for the update step, where ck is a constant and z1:k denotes the measurements from time instant 1 to k. The state transition prior p(xk∣xk-1) and the measurement likelihood p(zk∣xk) are obtained from the movement model in (1) and the measurement model in (2), respectively.


## 2.3. Channel-SLAM

In the following, we will revise the Channel-SLAM algorithm from [13, 17]. Figure 2 gives an overview of the two stages of Channel-SLAM. In the first stage, the parameters of the signal components received by the user via different propagation paths are estimated. The resulting estimates are used as measurement input in the second stage, where the states of the physical and virtual transmitters and the user position are estimated simultaneously in a SLAM scheme. Further sensors, such as an inertial measurement unit (IMU), may be included in the second stage. The locations of both the physical and virtual transmitters are assumed to be unknown. Thus, Channel-SLAM does not differentiate between physical and virtual transmitters, and the term transmitter comprises both physical and virtual transmitters in the following. Each signal component arriving at the receiver corresponds to one transmitter.Figure 2
Based on the received signal, the parameters of the propagation paths are estimated in the first step by the KEST algorithm. In the second step, the estimates serve as measurements for estimating the positions of the user and the physical and virtual transmitters. In addition, user heading rates of change measurements from an IMU are incorporated in the second step.The RF propagation channel between the physical transmitter and the user equipped with a receiver is assumed to be a linear and time-variant multipath channel.The received signal is modeled as a superposition of signal components of the transmit signals(t), where the jth signal component is defined by a complex amplitude aj(tk) and a delay dj(tk) at time tk. The signal received by the user at time instant tk is(5)yτ,tk=∑jajtksτ-djtk+nτ,where n(τ) is a sample from a colored noise sequence incorporating both dense multipath components (DMCs) and additive Gaussian noise. The channel is assumed to be constant during the short time interval when the received signal is sampled at time instant k.The physical transmitter continuously broadcasts the signals(t) that is known to the user. At the user side, the parameters of the signal components arriving at the receiver are estimated. Such parameters can in general be the complex amplitude, time of arrival (ToA), angle of arrival (AoA), or Doppler shift, depending on the available hardware and the scenario. For the signal parameter estimation, we use the KEST algorithm [36]. The KEST estimator works in two stages. In an inner stage, a ML parameter estimator such as Space-Alternating Generalized Expectation-Maximization (SAGE) [37] estimates the parameters of the signal components jointly on a snapshot basis. An outer stage tracks these estimated parameters over time with a Kalman filter and keeps track of the number of signal components. The KEST estimator is in general able to handle the DMCs in the noise term in (5). However, DMC handling is not implemented in our evaluations, leading to a model mismatch in KEST and hence to a higher variance in the parameter estimation. However, we do not expect many DMCs in our evaluation scenarios. In an indoor scenario, for example, DMCs need to be considered [38].In the literature, there are alternatives to the KEST estimator. For example, the authors of [39] track signal component parameters based on an EKF, though the authors of [40] showed that the KEST estimator is more robust in resolving signal components that are close to each other in the state space. In [41], an EKF is used as well for parameter estimation, while the position estimation is based on the time difference of arrival (TDoA) of virtual transmitters. The authors of [42] consider the linearization of the observation model in the EKF a major drawback that might lead to a tracking loss.In the second stage of Channel-SLAM, we use only the delays, that is, ToAs, and AoAs, estimated by KEST as measurement inputs. Hence, after sampling the received signal, the KEST estimates at time instantk are comprised in the vector(6)zk=dkTθkTT,where(7)dk=d1,k⋯dNTX,kTare the ToA estimates for the NTX signal components, or transmitters, and(8)θk=θ1,k⋯θNTX,kTare the corresponding AoA estimates. Note that the number of signal components and thus transmitters may change over time. Nevertheless, for notational convenience, we omit the time instant index k in NTX.In the second stage of Channel-SLAM, the user statexu,k is estimated simultaneously with the state of the transmitters xTX,k. The entire state vector is hence(9)xk=xu,kTxTX,kTT=xu,kTxTX,k1T⋯xTX,kNTXTT,where xTX,kj is the state of the jth transmitter. As we consider a two-dimensional scenario, the user state at time instant k is defined by(10)xu,k=xkykvx,kvy,kT=pu,kTvu,kTT,where the user position is defined by pu,k=[xkyk]T and the user velocity by vu,k=[vx,kvy,k]T. Each transmitter is defined by its location pTX,k=[xTX,kyTX,k]T and a clock offset τ0,k at time instant k. The state vector of the jth transmitter is hence defined by(11)xTX,kj=xTX,kjyTX,kjτ0,kjT=pTX,kjTτ0,kjT.Our goal is to find the minimum mean square error (MMSE) estimator forxk, which is defined as (12)x^k=∫xkpxk∣z1:kdxk, where z1:k denotes all measurements up to time instant k.We use a recursive Bayesian estimation scheme as in Section2.2 to estimate the posterior PDF p(xk∣z1:k). This posterior can be factorized as(13)pxk∣z1:k=pxTX,k,xu,k∣z1:k=pxTX,k∣xu,k,z1:kpxu,k∣z1:k.The signal components arriving at the receiver are assumed to be independent of each other; that is, we assume they interact with distinct objects. Assuming independence among the measurements for distinct transmitters (on the one hand, the parameters of the signal components are estimated jointly by the KEST algorithm, and hence these estimates might be correlated between signal components and between the parameters; on the other hand, the correlation is likely to have effect only on a short term basis as KEST estimates are unbiased when observed over a longer time), the first factor in (13) can be factorized further as(14)pxTX,0:k∣xu,0:k,z1:k=∏j=1NTXpxTX,0:kj∣xu,0:k,z1:kj.With the above factorization, the transmitter states are estimated independently from each other.As we consider a static scenario, the virtual transmitters are static as well, and the transition prior for thejth transmitter is calculated as(15)pxTX,kj∣xTX,k-1j=δxTX,kj-xTX,k-1j,where δ(·) denotes the Dirac distribution.For the prediction of the user, additional sensors such as an IMU carried by the user may be integrated into the movement model. Within this paper, we assume only heading change rate measurements from a gyroscope to be available, though, and no knowledge on the user speed. With the gyroscope heading change rateΔβ,k, we predict the movement of the user by(16)xu,k=12Tk1202RΔβ,kxu,k-1=Fu,kxu,k-1,where Tk denotes the time between instants k-1 and k. The two-dimensional rotation matrix R(Δβ,k) is defined as(17)RΔβ,k=cos⁡Δβ,k+wk-sin⁡Δβ,k+wksin⁡Δβ,k+wkcos⁡Δβ,k+wk,where wk is the heading noise which is distributed following a von Mises distribution. Hence, the function fk in (1) can be expressed in our case as(18)fkxk-1,vk-1=Fu,k04×3NTX03NTX×413NTXxk-1+vk-1,where the process noise covariance matrix Q is diagonal.As depicted in Figure3, an AoA measurement for a transmitter j describes the angle θj,k between the user heading direction vu,k and the incoming signal from the transmitter. The measurement noise for the ToA and AoA measurements is assumed to be zero-mean Gaussian distributed with variances σd,j2 and σθ,j2, respectively, for the jth transmitter. Also, we assume no cross-correlation between the single ToA and AoA measurements. The likelihood for the measurement vector zk conditioned on the state vector xk is therefore the product(19)pzk∣xk=∏j=1NTXNdj,k;d^j,k,σd,j2Nθj,k;θ^j,k,σθ,j2,where the predicted ToA between the user and the jth transmitter is(20)d^j,k=1c0pu,k-pTX,kj+τ0,kj,and the predicted AoA is calculated as(21)θ^j,k=atan2⁡yk-yTX,kj,xk-xTX,kj-atan2⁡vy,k,vx,k.The function atan2(y,x) calculates the four-quadrant inverse tangent function. It returns the counterclockwise angle between the positive x-axis and the point given by the coordinates (x,y).Figure 3
The user moves in the directionvu,k at time instant k. The heading change rate from the IMU is Δβ,k. The ToA and AoA measurements for the signal from the jth transmitter are dj,k and θj,k, respectively, where θj,k describes the angle between the heading direction vu,k of the user and the arriving signal.


## 3. Nonlinear Kalman Filtering

### 3.1. Unscented Transform

If a random variablex is transformed by a function g(·) such that y=g(x), the statistics of y cannot always be calculated in closed form. Monte Carlo (MC) methods try to estimate the statistics of y from a set of randomly chosen sample points of x that undergo the transformation g(·). For the unscented transform, a set of the so-called sigma points is propagated through the function g(·) to obtain transformed sigma points yielding the statistics of y. However, the sigma points are not chosen randomly, but in a deterministic manner, which is the fundamental difference to MC methods.Based on the unscented transform, numerical approximations of integrals can be derived. In particular, for the case when the integrand is a product of an arbitrary functiong(x) of the integration variable x and a Gaussian PDF Nx;μx,Cx an integration rule of the form(22)∫gxNx;μx,Cxdx≈∑m=1NsigωmgXmcan be applied, where Xm is the mth of the Nsig sigma points with its associated weight ωm. The idea of the UKF is to approximate the posterior PDF in recursive Bayesian estimation by a Gaussian PDF. Hence, the integral in the prediction step is approximated by the integration rule in (22). The authors of [43] provide further insight into sigma point methods and their relation to Gaussian process quadrature.


### 3.2. Choice of Sigma Points

In the literature, different sets of sigma points have been proposed for the unscented transform [44]. Let Xm be the mth sigma point and ωm its associated weight. The dimension, mean, and covariance of the random variable x are denoted by N, μx, and Cx, respectively. In [23], the sigma points and their weights are defined for some κ∈R as(23)X0=μx,ω0=κκ+N,Xm=μx+N+κCxm,ωm=12κ+N,Xm+N=μx-N+κCxm,ωm+N=12κ+N,where m=1,…,N, (A)m denotes the mth row or column of the matrix A, and (N+κ)Cx is factorized into(24)N+κCx=N+κCxN+κCxT.This definition leads toNsig=2N+1 sigma points. The authors of [45] presented the cubature Kalman filter (CKF) with an intuitive derivation of the choice of sigma points and their weights. The CKF differs from the UKF only in the choice of the sigma points. Its derivation is based on the fact that the approximation of an integral using the unscented transform as in (22) is exact for g(x) being a monomial of an order not greater than some integer d. The resulting sigma points are the points in (23) for κ=0. Since the weight of the first sigma point is zero, there are only 2N effective points. Although the derivation in [45] gives useful insight into the UKF, the same choice of sigma points had been proposed in [44] before.The Appendix summarizes the equations for the prediction and update steps of the UKF. If the state transition model in (1) or the measurement model in (2) are linear or if Gaussian noise is assumed in the state transition or measurement model, methods from [46] can be applied to decrease the computational complexity of the UKF.


## 3.1. Unscented Transform

If a random variablex is transformed by a function g(·) such that y=g(x), the statistics of y cannot always be calculated in closed form. Monte Carlo (MC) methods try to estimate the statistics of y from a set of randomly chosen sample points of x that undergo the transformation g(·). For the unscented transform, a set of the so-called sigma points is propagated through the function g(·) to obtain transformed sigma points yielding the statistics of y. However, the sigma points are not chosen randomly, but in a deterministic manner, which is the fundamental difference to MC methods.Based on the unscented transform, numerical approximations of integrals can be derived. In particular, for the case when the integrand is a product of an arbitrary functiong(x) of the integration variable x and a Gaussian PDF Nx;μx,Cx an integration rule of the form(22)∫gxNx;μx,Cxdx≈∑m=1NsigωmgXmcan be applied, where Xm is the mth of the Nsig sigma points with its associated weight ωm. The idea of the UKF is to approximate the posterior PDF in recursive Bayesian estimation by a Gaussian PDF. Hence, the integral in the prediction step is approximated by the integration rule in (22). The authors of [43] provide further insight into sigma point methods and their relation to Gaussian process quadrature.


## 3.2. Choice of Sigma Points

In the literature, different sets of sigma points have been proposed for the unscented transform [44]. Let Xm be the mth sigma point and ωm its associated weight. The dimension, mean, and covariance of the random variable x are denoted by N, μx, and Cx, respectively. In [23], the sigma points and their weights are defined for some κ∈R as(23)X0=μx,ω0=κκ+N,Xm=μx+N+κCxm,ωm=12κ+N,Xm+N=μx-N+κCxm,ωm+N=12κ+N,where m=1,…,N, (A)m denotes the mth row or column of the matrix A, and (N+κ)Cx is factorized into(24)N+κCx=N+κCxN+κCxT.This definition leads toNsig=2N+1 sigma points. The authors of [45] presented the cubature Kalman filter (CKF) with an intuitive derivation of the choice of sigma points and their weights. The CKF differs from the UKF only in the choice of the sigma points. Its derivation is based on the fact that the approximation of an integral using the unscented transform as in (22) is exact for g(x) being a monomial of an order not greater than some integer d. The resulting sigma points are the points in (23) for κ=0. Since the weight of the first sigma point is zero, there are only 2N effective points. Although the derivation in [45] gives useful insight into the UKF, the same choice of sigma points had been proposed in [44] before.The Appendix summarizes the equations for the prediction and update steps of the UKF. If the state transition model in (1) or the measurement model in (2) are linear or if Gaussian noise is assumed in the state transition or measurement model, methods from [46] can be applied to decrease the computational complexity of the UKF.


## 4. Derivation of the Gaussian Sum Particle Filter

### 4.1. The Rao-Blackwellized Gaussian Sum Particle Filter

The factorization in (13) allows for estimating the user state independently from the transmitter states. For the estimation of the user state in the RBGSPF, we use a SIR particle filter [26, 47]. The single transmitter states xTX,kj are estimated independently from each other following (14). Each transmitter state is represented by a Gaussian mixture model or Gaussian sum model [48]. The posterior PDF of each of the NUKF Gaussian components in a Gaussian mixture is estimated by a UKF. The structure of the resulting RBGSPF representation is shown in Figure 4.Figure 4
Structure of the RBGSPF representation: the user statexu is represented by a number of particles. Each particle estimates the transmitters’ states on its own. Each of the NTX transmitter states is represented by a sum of NUKF Gaussian distributions Nl with associated weights wl.A particle filter is a MC based method, where the posterior PDF is represented by a number of random samples, called particles, with associated weights. The user posterior PDF is approximated as(25)pxu,k∣z1:k=∑i=1Npwkiδxu,k-xu,ki,where xu,ki is the ith user particle, wki its associated weight, and Np the number of particles in the particle filter. From the structure of (14), we see that the transmitter states are estimated for each user particle independently.The posterior distribution of the state of each transmitter is approximated by a Gaussian mixture model. In a Gaussian mixture model, a PDF is described as a sum of weighted Gaussian PDFs, each described by a mean and a covariance. Hence, the posterior PDF of the state of thejth transmitter of the ith user particle is represented as [48](26)pxTX,ki,j∣xu,ki,z1:kj=∑l=1NUKFwki,j,l×NxTX,ki,j;xTX,k∣ki,j,l,Pk∣ki,j,l,where z1:kj are the measurements for the jth transmitter and xTX,k∣ki,j,l, Pk∣ki,j,l, and wki,j,l are the mean, the covariance matrix, and the weight, respectively, of the lth Gaussian component of the Gaussian mixture for the jth transmitter. Both xTX,k∣ki,j,l and Pk∣ki,j,l are obtained from the update step of the corresponding UKF. Similarly, the likelihood for the measurement of the jth transmitter of the ith user particle is(27)pzkj∣xu,ki,xTX,ki,j=∑l=1NUKFwki,j,l×Nzkj;z^ki,j,l,Rkj,where Rkj is the measurement noise covariance matrix for the jth transmitter. The number NUKF of Gaussian components might differ between transmitters, user particles, and time instants. However, for notational convenience, we drop the particle, transmitter, and time instant indices of NUKF. The predicted measurement of the lth Gaussian component for the jth transmitter of the ith user particle,(28)z^ki,j,l=d^ki,j,lθ^ki,j,lT,consists of the predicted ToA measurement(29)d^ki,j,l=1c0pu,ki-pTX,ki,j,l+τ0,ki,j,land the predicted AoA measurement(30)θ^ki,j,l=atan2⁡yki-yTX,ki,j,l,xki-xTX,ki,j,l-atan2⁡vy,ki,vx,ki,where the ith user particle is (31)xu,ki=xkiykivx,kivy,kiT=pu,kiTvu,kiTT, and the mean of the corresponding lth Gaussian component of the jth transmitter is (32)xTX,k∣ki,j,l=xTX,ki,j,lyTX,ki,j,lτ0,ki,j,lT=pTX,ki,j,lTτ0,ki,j,lT.In the prediction step of the user particle filter, new particles are sampled based on the transition priorp(xu,k∣xu,k-1). Hence, the particle xu,ki is drawn as(33)xu,ki=fu,kxu,k-1i,vu,k-1,where the function fu,k· describes the user movement model, and vu,k-1 is a noise sample drawn from the user process noise PDF. For the prediction and the update step of a Gaussian component of a transmitter’s state, the equations of the UKF are summarized in the Appendix.


### 4.2. Derivation of the Particle Weight Calculation

In the following, we will derive the calculation of the particle weights in the user particle filter and of the weights for the Gaussian components in the Gaussian mixture models used to describe the PDFs of the transmitter states. As the importance density of the SIR particle filter is the state transition prior and resampling of the particles is performed at every time instant, the weight for theith particle at time instant k is given by [26] (34)wki∝pzk∣xu,0:ki,z1:k-1. This expression can be written as(35)wki∝∫pzk∣xu,0:ki,xTX,ki,z1:k-1×pxTX,ki∣xu,0:ki,z1:k-1dxTX,ki∝∫pzk∣xu,ki,xTX,ki×pxTX,ki∣xu,ki,zk-1dxTX,ki,where we use the assumption of a first-order hidden Markov model. With the assumption that the measurements are independent for different transmitters, (35) can be expressed as(36)wki∝∏j=1NTX∫pzkj∣xu,ki,xTX,ki,j×pxTX,ki,j∣xu,ki,zk-1jdxTX,ki,j.Furthermore, using the Gaussian mixture model representation from (26) and (27) and assuming Gaussian measurement noise, the integrand can be expressed as a sum of weighted Gaussian PDFs; namely,(37)pzkj∣xu,ki,xTX,ki,jpxTX,ki,j∣xu,ki,zk-1j=∑l=1NUKFwk-1i,j,lpz,ki,j,lpTX,ki,j,l.In (37), we defined for notational brevity(38)pz,ki,j,l=Nzkj;z^ki,j,l,Rkj,(39)pTX,ki,j,l=NxTX,ki,j,l;xTX,k∣ki,j,l,Pk∣ki,j,l,where Rkj is the measurement noise covariance matrix for the jth transmitter and xTX,ki,j,l denotes the state of the lth Gaussian component of the jth transmitter of the ith user particle. As we assume no correlation among the ToA and AoA measurements, we have(40)Rkj=σd,j200σθ,j2.Inserting (37) into (36) leads to(41)wki∝∏j=1NTX∑l=1NUKFwk-1i,j,l∫pz,ki,j,lpTX,ki,j,ldxTX,ki,j,l.The predicted measurementz^ki,j,l in (38), defined in (28), is a nonlinear function of xTX,ki,j,l. We express the measurement likelihood in (38) explicitly as a function g(·) of xTX,ki,j,l, resulting in(42)gxTX,ki,j,l=pz,ki,j,l=Ndj,k;d^ki,j,l,σd,j2×Nθj,k;θ^ki,j,l,σθ,j2.Due to the nonlinearity in (42), stemming from (29) and (30), the integral in (41) cannot be solved analytically. Instead, we use the approximation from (22) to calculate the integral as(43)∫gxTX,ki,j,lNxTX,ki,j,l;xTX,k∣ki,j,l,Pk∣ki,j,ldxTX,ki,j,l≈∑m=1NsigωmgXm,where the sigma points Xm and their weights ωm can be calculated by (23), where (44)μx=xTX,k∣ki,j,l,Cx=Pk∣ki,j,l, and N is the dimension of a transmitter’s state; that is, N=3. Finally, the weight of the ith particle is calculated as(45)wki∝∏j=1NTX∑l=1NUKFwk-1i,j,l∑m=1NsigωmgXm.It follows directly from (45) that the weights of the Gaussian components are updated by(46)wki,j,l∝wk-1i,j,l∑m=1NsigωmgXm.Note that the weightswki and wki,j,l in (45) and (46) of the user particles and the Gaussian components, respectively, are not yet normalized. Since resampling is performed at every time instant in the SIR particle filter, the user particle weights wki do not depend on the weights wk-1i from the previous time instant [26].


### 4.3. Merging and Pruning of Gaussian Components

When the KEST estimator detects a new signal component, a new transmitter is initialized for each user particle based on the ToA and AoA measurement for that new transmitter at the current time instant. The posterior PDF of the new transmitter is represented by a number of Gaussian PDFs, whose means are initially placed on a grid dependent on the measurement. The number of Gaussian components depends on the measurement as well.As the user travels through a scenario, the means, covariances, and weights of the Gaussian components of a transmitter’s state posterior PDF change over time depending on the available measurements. The mean and covariance of a Gaussian component may be regarded as a hypothesis and a corresponding uncertainty, respectively, for the state of a transmitter. If the weight of a Gaussian component becomes smaller, the hypothesis for that state of the transmitter becomes less likely. Hence, if the weight of a Gaussian component falls below a thresholdρ, the Gaussian component is pruned; that is, its weight wk-1i,j,l is set to zero. If the means of two Gaussian components get very close to each other, they may be merged in order to reduce the computational complexity.The final algorithm for one time instantk>0 of the RBGSPF is summarized in Algorithm 1. For a particle filter resampling algorithm, we refer to [26]. Note again that we have dropped the particle and transmitter indices in NUKF.Algorithm 1:RBGSPF for time instant k>0.
Input:  xu,k-1i, xTX,k-1i,j,l and wk-1i,j,l for i=1,…,Np, j=1,…,NTX, l=1,…,NUKF, zk Output:  xu,ki, xTX,ki,j,l and wki,j,l for i=1,…,Np, j=1,…,NTX, l=1,…,NUKF (1) fori=1,…,Npdo (2)draw new user particle xu,ki using (33); (3)ifany new signal components detectedthen (4)initialize the new transmitter(s) based on zk; (5)iftrack of any signal components lostthen (6)delete the corresponding transmitter(s); (7)forj=1,…,NTXdo (8)forl=1,…,NUKFdo (9)perform UKF prediction and update to calculate xTX,ki,j,l using the UKF equations in the Appendix; (10)calculate the weight wki,j,l with (46); (11)ifwki,j,l<ρthen (12)prune this Gaussian component: set wki,j,l=0; (13)calculate the weight wki with (45); (14) fori=1,…,Npdo (15)forj=1,…,NTXdo (16)normalize the weights wki,j,l for l=1,…,NUKF; (17) normalize the weights wki for i=1,…,Np; (18) resample the user particles xu,ki [26];


## 4.1. The Rao-Blackwellized Gaussian Sum Particle Filter

The factorization in (13) allows for estimating the user state independently from the transmitter states. For the estimation of the user state in the RBGSPF, we use a SIR particle filter [26, 47]. The single transmitter states xTX,kj are estimated independently from each other following (14). Each transmitter state is represented by a Gaussian mixture model or Gaussian sum model [48]. The posterior PDF of each of the NUKF Gaussian components in a Gaussian mixture is estimated by a UKF. The structure of the resulting RBGSPF representation is shown in Figure 4.Figure 4
Structure of the RBGSPF representation: the user statexu is represented by a number of particles. Each particle estimates the transmitters’ states on its own. Each of the NTX transmitter states is represented by a sum of NUKF Gaussian distributions Nl with associated weights wl.A particle filter is a MC based method, where the posterior PDF is represented by a number of random samples, called particles, with associated weights. The user posterior PDF is approximated as(25)pxu,k∣z1:k=∑i=1Npwkiδxu,k-xu,ki,where xu,ki is the ith user particle, wki its associated weight, and Np the number of particles in the particle filter. From the structure of (14), we see that the transmitter states are estimated for each user particle independently.The posterior distribution of the state of each transmitter is approximated by a Gaussian mixture model. In a Gaussian mixture model, a PDF is described as a sum of weighted Gaussian PDFs, each described by a mean and a covariance. Hence, the posterior PDF of the state of thejth transmitter of the ith user particle is represented as [48](26)pxTX,ki,j∣xu,ki,z1:kj=∑l=1NUKFwki,j,l×NxTX,ki,j;xTX,k∣ki,j,l,Pk∣ki,j,l,where z1:kj are the measurements for the jth transmitter and xTX,k∣ki,j,l, Pk∣ki,j,l, and wki,j,l are the mean, the covariance matrix, and the weight, respectively, of the lth Gaussian component of the Gaussian mixture for the jth transmitter. Both xTX,k∣ki,j,l and Pk∣ki,j,l are obtained from the update step of the corresponding UKF. Similarly, the likelihood for the measurement of the jth transmitter of the ith user particle is(27)pzkj∣xu,ki,xTX,ki,j=∑l=1NUKFwki,j,l×Nzkj;z^ki,j,l,Rkj,where Rkj is the measurement noise covariance matrix for the jth transmitter. The number NUKF of Gaussian components might differ between transmitters, user particles, and time instants. However, for notational convenience, we drop the particle, transmitter, and time instant indices of NUKF. The predicted measurement of the lth Gaussian component for the jth transmitter of the ith user particle,(28)z^ki,j,l=d^ki,j,lθ^ki,j,lT,consists of the predicted ToA measurement(29)d^ki,j,l=1c0pu,ki-pTX,ki,j,l+τ0,ki,j,land the predicted AoA measurement(30)θ^ki,j,l=atan2⁡yki-yTX,ki,j,l,xki-xTX,ki,j,l-atan2⁡vy,ki,vx,ki,where the ith user particle is (31)xu,ki=xkiykivx,kivy,kiT=pu,kiTvu,kiTT, and the mean of the corresponding lth Gaussian component of the jth transmitter is (32)xTX,k∣ki,j,l=xTX,ki,j,lyTX,ki,j,lτ0,ki,j,lT=pTX,ki,j,lTτ0,ki,j,lT.In the prediction step of the user particle filter, new particles are sampled based on the transition priorp(xu,k∣xu,k-1). Hence, the particle xu,ki is drawn as(33)xu,ki=fu,kxu,k-1i,vu,k-1,where the function fu,k· describes the user movement model, and vu,k-1 is a noise sample drawn from the user process noise PDF. For the prediction and the update step of a Gaussian component of a transmitter’s state, the equations of the UKF are summarized in the Appendix.


## 4.2. Derivation of the Particle Weight Calculation

In the following, we will derive the calculation of the particle weights in the user particle filter and of the weights for the Gaussian components in the Gaussian mixture models used to describe the PDFs of the transmitter states. As the importance density of the SIR particle filter is the state transition prior and resampling of the particles is performed at every time instant, the weight for theith particle at time instant k is given by [26] (34)wki∝pzk∣xu,0:ki,z1:k-1. This expression can be written as(35)wki∝∫pzk∣xu,0:ki,xTX,ki,z1:k-1×pxTX,ki∣xu,0:ki,z1:k-1dxTX,ki∝∫pzk∣xu,ki,xTX,ki×pxTX,ki∣xu,ki,zk-1dxTX,ki,where we use the assumption of a first-order hidden Markov model. With the assumption that the measurements are independent for different transmitters, (35) can be expressed as(36)wki∝∏j=1NTX∫pzkj∣xu,ki,xTX,ki,j×pxTX,ki,j∣xu,ki,zk-1jdxTX,ki,j.Furthermore, using the Gaussian mixture model representation from (26) and (27) and assuming Gaussian measurement noise, the integrand can be expressed as a sum of weighted Gaussian PDFs; namely,(37)pzkj∣xu,ki,xTX,ki,jpxTX,ki,j∣xu,ki,zk-1j=∑l=1NUKFwk-1i,j,lpz,ki,j,lpTX,ki,j,l.In (37), we defined for notational brevity(38)pz,ki,j,l=Nzkj;z^ki,j,l,Rkj,(39)pTX,ki,j,l=NxTX,ki,j,l;xTX,k∣ki,j,l,Pk∣ki,j,l,where Rkj is the measurement noise covariance matrix for the jth transmitter and xTX,ki,j,l denotes the state of the lth Gaussian component of the jth transmitter of the ith user particle. As we assume no correlation among the ToA and AoA measurements, we have(40)Rkj=σd,j200σθ,j2.Inserting (37) into (36) leads to(41)wki∝∏j=1NTX∑l=1NUKFwk-1i,j,l∫pz,ki,j,lpTX,ki,j,ldxTX,ki,j,l.The predicted measurementz^ki,j,l in (38), defined in (28), is a nonlinear function of xTX,ki,j,l. We express the measurement likelihood in (38) explicitly as a function g(·) of xTX,ki,j,l, resulting in(42)gxTX,ki,j,l=pz,ki,j,l=Ndj,k;d^ki,j,l,σd,j2×Nθj,k;θ^ki,j,l,σθ,j2.Due to the nonlinearity in (42), stemming from (29) and (30), the integral in (41) cannot be solved analytically. Instead, we use the approximation from (22) to calculate the integral as(43)∫gxTX,ki,j,lNxTX,ki,j,l;xTX,k∣ki,j,l,Pk∣ki,j,ldxTX,ki,j,l≈∑m=1NsigωmgXm,where the sigma points Xm and their weights ωm can be calculated by (23), where (44)μx=xTX,k∣ki,j,l,Cx=Pk∣ki,j,l, and N is the dimension of a transmitter’s state; that is, N=3. Finally, the weight of the ith particle is calculated as(45)wki∝∏j=1NTX∑l=1NUKFwk-1i,j,l∑m=1NsigωmgXm.It follows directly from (45) that the weights of the Gaussian components are updated by(46)wki,j,l∝wk-1i,j,l∑m=1NsigωmgXm.Note that the weightswki and wki,j,l in (45) and (46) of the user particles and the Gaussian components, respectively, are not yet normalized. Since resampling is performed at every time instant in the SIR particle filter, the user particle weights wki do not depend on the weights wk-1i from the previous time instant [26].


## 4.3. Merging and Pruning of Gaussian Components

When the KEST estimator detects a new signal component, a new transmitter is initialized for each user particle based on the ToA and AoA measurement for that new transmitter at the current time instant. The posterior PDF of the new transmitter is represented by a number of Gaussian PDFs, whose means are initially placed on a grid dependent on the measurement. The number of Gaussian components depends on the measurement as well.As the user travels through a scenario, the means, covariances, and weights of the Gaussian components of a transmitter’s state posterior PDF change over time depending on the available measurements. The mean and covariance of a Gaussian component may be regarded as a hypothesis and a corresponding uncertainty, respectively, for the state of a transmitter. If the weight of a Gaussian component becomes smaller, the hypothesis for that state of the transmitter becomes less likely. Hence, if the weight of a Gaussian component falls below a thresholdρ, the Gaussian component is pruned; that is, its weight wk-1i,j,l is set to zero. If the means of two Gaussian components get very close to each other, they may be merged in order to reduce the computational complexity.The final algorithm for one time instantk>0 of the RBGSPF is summarized in Algorithm 1. For a particle filter resampling algorithm, we refer to [26]. Note again that we have dropped the particle and transmitter indices in NUKF.Algorithm 1:RBGSPF for time instant k>0.
Input:  xu,k-1i, xTX,k-1i,j,l and wk-1i,j,l for i=1,…,Np, j=1,…,NTX, l=1,…,NUKF, zk Output:  xu,ki, xTX,ki,j,l and wki,j,l for i=1,…,Np, j=1,…,NTX, l=1,…,NUKF (1) fori=1,…,Npdo (2)draw new user particle xu,ki using (33); (3)ifany new signal components detectedthen (4)initialize the new transmitter(s) based on zk; (5)iftrack of any signal components lostthen (6)delete the corresponding transmitter(s); (7)forj=1,…,NTXdo (8)forl=1,…,NUKFdo (9)perform UKF prediction and update to calculate xTX,ki,j,l using the UKF equations in the Appendix; (10)calculate the weight wki,j,l with (46); (11)ifwki,j,l<ρthen (12)prune this Gaussian component: set wki,j,l=0; (13)calculate the weight wki with (45); (14) fori=1,…,Npdo (15)forj=1,…,NTXdo (16)normalize the weights wki,j,l for l=1,…,NUKF; (17) normalize the weights wki for i=1,…,Np; (18) resample the user particles xu,ki [26];


## 5. Data Association

Data association is of crucial importance for robust long-term SLAM. It describes the correspondences among landmarks, which are transmitters in multipath assisted positioning. In Figure5, a user travels along its trajectory. In Region (I), the LoS signal from the transmitter Tx can be tracked. This signal is lost in Region (II) and regained in Region (III). However, KEST is not able to retrack a former path. Hence, when the user enters Region (III), KEST detects a new signal component, and consequently a new transmitter is initialized. However, the transmitter is the same as that which had been observed in Region (I).Figure 5
A user moves along the trajectory. The LoS signal to the transmitter in Region (I) is lost in Region (II) temporarily due to blocking by an obstacle and received again in Region (III).We define the set of transmitters that had been observed previously but are not detected any more, asoldtransmitters. When a new signal component is detected by KEST, a new transmitter has to be initialized. Consequently, two cases may arise:(1)
the new signal component corresponds indeed to a new transmitter, or(2)
the new signal component corresponds to an old transmitter that had been observed before already.Data association is the decision on the above two cases when a new signal component is detected. In the first case, a new transmitter is initialized for the newly detected signal component. In the second case, the newly detected signal component is associated with a previously observable, that is, old, transmitter.In [33], a multiple hypothesis tracking (MHT) association method was introduced and derived for FastSLAM, where the user state is represented by a particle filter and each landmark state by an EKF. In [19], the same method has been derived for a Rao-Blackwellized particle filter. In the following, we will derive the method for the RBGSPF derived in Section 4.Each user particle decides for associations individually and thus carries a hypothesis for associations. Hence, association decisions are hard decisions for each particle. Regarding the ensemble of user particles, though, there are many different hypotheses on associations in the user state estimate, and the association method can be regarded as a soft decision method. Consequently, the state vector of the user is increased by data association.In the following, we describe how to make an association decision for a single user particle, where we omit the particle indexi in the association variables for notational brevity. The value of the association variable nk denotes an association of the new transmitter with the old transmitter nk. We denote the marginalized likelihood of the measurement of the new transmitter that is to be initialized at time instant k by pnk assuming that the new transmitter is associated with the old transmitter nk. The set of association decisions up to time instant k-1 is denoted by Nk-1. From [33], we have(47)pnk=pzk∣nk,Nk-1,xu,ki,z1:k-1=∫pzk∣xTX,ki,nk,nk,Nk-1,xu,ki,z1:k-1×pxTX,ki,nk∣nk,Nk-1,xu,ki,z1:k-1dxTX,ki,nk,where xTX,ki,nk denotes the state vector of the nkth transmitter for the ith user particle.Assuming a first-order hidden Markov model, the first integrand in (47) can be simplified to(48)pzk∣xTX,ki,nk,nk,Nk-1,xu,ki,z1:k-1=pzk∣xTX,ki,nk,nk,Nk-1,xu,ki.Since we use a Gaussian mixture model to represent the single transmitter states, the second integrand in (47) can be rewritten as(49)pxTX,ki,nk∣nk,Nk-1,xu,ki,z1:k-1=∑l=1NUKFwki,nk,lNxTX,ki,nk;xTX,k∣k-1i,nk,l,Pk∣k-1i,nk,l.Inserting (48) and (49) into (47) yields(50)pnk=∑l=1NUKFwki,nk,l∫pz,ki,nk,lpTX,ki,nk,ldxTX,ki,nk,l,where pz,ki,nk,l is defined as in (38), and(51)pTX,k∣k-1i,nk,l=NxTX,ki,nk,l;xTX,k∣k-1i,nk,l,Pk∣k-1i,nk,l.Similar to (42), we define(52)gxTX,ki,nk,l=pz,ki,nk,l=Ndnk,k;d^ki,nk,l,σd,nk2×Nθnk,k;θ^ki,nk,l,σθ,nk2, which is nonlinear in xTX,ki,nk,l, and rewrite the integral in (50) as (53)∫pz,ki,nk,lpTX,k∣k-1i,nk,ldxTX,ki,nk,l=∫gxTX,ki,nk,l×NxTX,ki,nk,l;xTX,k∣k-1i,nk,l,Pk∣k-1i,nk,ldxTX,ki,nk,l.Approximating the integral using (22) and inserting it into (50) finally yield(54)pnk=∑l=1NUKFwki,nk,l∑m=1NsigωmgXm.The sigma points Xm and their weights ωm can again be calculated by (23) with (55)μx=xTX,k∣k-1i,nk,l,Cx=Pk∣k-1i,nk,l.The authors of [33] propose two ways to come to an association decision, a ML method and data association sampling (DAS). The probability for making no association is defined and denoted by p0. The set of indices of old transmitters that have not yet been and hence might be associated is denoted by Γk.For the ML association method, the association of the new transmitter with the old transmitternk is chosen to be(56)n^ML,k=argmaxnk∈Γk∪0⁡pnk.In DAS, an association is sampled based on the likelihoods cpnk for nk∈Γk∪0, where c is a normalization constant.If a detected transmitter is associated with an old transmitter, the new transmitter can be initialized with the posterior PDF of the associated old transmitter. Thus, data association has to be incorporated in Line (4) of Algorithm1.The above method describes how to take association decisions if no more than one new transmitter is initialized at one time instant, that is, if no more than one new signal component is detected by KEST at a time instant. In case of multiple transmitters being initialized at the same time instant, a greedy algorithm [20] may be applied.


## 6. Evaluations

In the following, we evaluate the RBGSPF derived in Section4 by means of simulations and actual outdoor measurements. For the evaluations, we implemented a square-root version of the cubature Kalman filter as in [45] for numerical stability. The sigma points are the ones in (23) for κ=0. Since the movement model of the transmitters is linear and we assume Gaussian noise, the prediction step can be calculated analytically. For the description of a prediction step of a square-root version of the conventional Kalman filter, we refer to [48].


### 6.1. Simulations in an Urban Scenario

A top view of the urban simulation scenario is depicted in Figure6. The thick black lines represent walls, for example, from buildings, that reflect RF signals, and the black circles are objects such as traffic light poles acting as scatterers. There is one physical transmitter in the scenario marked by the red upward triangle labeled Tx. The user travels with a constant speed of 10 m/s along the blue line with a loop around the central building. The initial and final user positions are labeled START and END, respectively. The traveled distances of the user are marked for every 50 m.Figure 6
The simulation scenario with thick black lines representing reflecting walls and black circles representing objects that scatter the RF signals. The physical transmitter is marked by the red upward triangle labeled Tx. The user travels along the blue line from START to END with one loop around the central building.The transmitter continuously broadcasts a signal that is known to the user and has a rectangular shape in frequency domain with a center frequency of 1.5 GHz and a bandwidth of 100 MHz. As we know the environment, a CIR and the received signal can be modeled for every user position with a simple ray-tracing approach. We incorporate first- and second-order reflections and scattering, that is, single and double reflections and/or scattering. The power loss for the signal being reflected is 3 dB and 6 dB when the signal is scattered at a point scatterer. The average signal-to-noise ratio (SNR) at the user is 7 dB.The user is equipped with an RF receiver and a two-dimensional, rectangular antenna array consisting of nine elements. Hence, both the ToA and the AoA estimates from KEST are incorporated in the estimation of the user and the transmitters’ states. Based on the received signal, KEST estimates the ToAs and AoAs every 50 ms.The results of the KEST estimator are plotted in Figure7. It shows the propagation distance, which is the ToA multiplied by the speed of light, of the signal components versus the traveled distance of the user. Each continuous line represents one signal component and its evolution as the user travels through the scenario. The color of each line indicates the normalized absolute value of the amplitude of the corresponding signal component in linear domain. Since signal components that are observable for a long time can contribute much better to Channel-SLAM than components which are observable only for a short time, only signal components that are observable for a user traveled distance of at least 35 m are plotted and used. Using all detected signal components would dramatically increase the computational complexity and hardly increase the positioning performance.Figure 7
The results of the KEST estimator for the simulations showing the propagation distances of signal components versus the user traveled distance. The propagation distances are the ToA multiplied by the speed of light. Only signal components that are observable for a traveled distance of at least 35 m are shown. The color indicates the normalized amplitude in linear domain.In Channel-SLAM, the user position is estimated relative to the physical and virtual transmitters in the scenario. Thus, to create a local coordinate system, the initial state of the user is assumed to be known. However, no prior knowledge on any transmitter is assumed. In the Rao-Blackwellized particle filter, the number of user particles in the particle filter is 4000, while the number of Gaussian components for each transmitter depends on the first ToA measurement for that transmitter.The root mean square error (RMSE) of the user position versus its traveled distance is plotted in Figure8. The red curve shows the RMSE if no associations among transmitters are made; that is, every signal component that is detected by the KEST algorithm is assumed to be a new transmitter. The RMSEs with the ML association method and DAS from Section 5 being applied are plotted in blue and green, respectively. Since the particle filter is a MC based method, all RMSE curves are averaged over 100 simulations.Figure 8
The RMSE of the user position versus the user traveled distance for the simulations. The red curve shows the RMSE if no associations among transmitters are made, the blue curve if the ML method for associations is applied, and the green curve for using DAS.As we assume the starting position of the user to be known, all three curves start with a low RMSE that increases linearly during the first 200 m as expected. The increase of the RMSE is less due to a bias in the position estimate but more due to an increasing uncertainty, that is, variance, about the user position. After approximately 200 m, the RMSE tends to decrease for all three curves. As more and more transmitters are observed, the weight for some user particles becomes small, and these particles are unlikely to be resampled. Towards the end of the track, the geometrical delusion of precision (GDOP) causes an increase in the RMSE, since most of the transmitters are observed from the same direction. After a traveled distance of around 370 m, several transmitters that had been observed in the beginning are observed again, and correspondences among them can be found. If data association methods are used, the RMSE decreases particularly in that region. The ML method and DAS show a similar performance. Note that there are several reasons for which associations among transmitters can be found. Examples are signal blocking or the geometry of the environment causing virtual transmitters to be observable only from certain regions. In addition, when KEST loses and regains track of a signal component if its received power fluctuates or if another signal component arrives at the receiver with a very small difference in delay, transmitters may be discarded and initialized again at a later point, and associations among them may be found. This explains the increasing positioning performance gain after approx. 50 m using data association.


### 6.2. Outdoor Measurements

In addition to the simulations as described above, we performed outdoor measurements on an airfield. A top view of the measurement scenario is depicted in Figure9. The grey area is an airplane hangar with solid metallic doors. The user track with a total length of 112.5 m is plotted in blue. The user walked along the track starting from the light blue cross labeled START to the black cross labeled END. The traveled distance of the user is marked after 25 m, 75 m, and 100 m. There is one physical transmitter marked by the red upward triangle labeled Tx. The user is in LoS to the physical transmitter throughout the entire track.Figure 9
Top view of the measurement scenario in front of a hangar. The physical transmitter location is marked by the red triangle labeled Tx. The user travels along the blue track from START to END. The expected virtual transmitter locations are marked by the magenta downward triangles.In the scenario, we have three fences labeled Fence 1, Fence 2, and Fence 3. We expect these fences and the hangar door to reflect the RF signal emitted by the physical transmitter. Hence, we expect a virtual transmitter for each of the fences and for the hangar door following Section2.1. The virtual transmitter corresponding to the reflection of the signal at Fence 1 is the magenta downward triangle labeled vTx2. It is located at the physical transmitter position mirrored at Fence 1. Likewise, the location of the virtual transmitter corresponding to Fence 2 is labeled as vTx3. For the reflection of the signal at the hangar doors, we expect the virtual transmitter located at the magenta triangle labeled as vTx1. The expected virtual transmitter corresponding to Fence 3, vTx4, is outside of the boundaries of Figure 9.The Medav RUSK broadband channel sounder [49] was used to perform the measurements. The transmit signal is a multitone signal with a center frequency of 1.51 GHz and a bandwidth of 100 MHz. The signal has 1281 subcarriers with equal gains and a total transmit power of 10 mW.The user was equipped with an RF receiver, recording a snapshot of the received signal every 1.024 ms. For later evaluation, the user carried a prism mounted next to the receiver antenna that was tracked by a tachymeter (Leica Geosystems TCRP1200) to obtain the ground truth of the user location in centimeter accuracy. In addition, the user carried an Xsens MTI-G-700 IMU. Only heading change rate measurements were used from the IMU.On both transmitter and receiver side, single antennas were used. Hence, no AoA information about the impinging signal components can be used for Channel-SLAM. Instead, only ToA estimates from KEST are incorporated. The likelihood function in (19) is adapted accordingly for the evaluation.The results of the KEST estimator for the outdoor measurements are plotted in Figure10. The colors indicate the power estimated by KEST in dBm. As for the simulations, only signal components that are observable for a long user traveled distance are plotted and used. In addition, the ground truth geometrical line-of-sight (GLoS) propagation distances from the physical and the expected virtual transmitters as in Figure 9 to the user are plotted by black lines. They match the KEST estimates very well, justifying the signal model in (5) without considering DMCs in KEST for the measurement scenario.Figure 10
The results of the KEST estimator for the outdoor measurements showing the propagation distances, that is, the ToAs multiplied by the speed of light, of signal components, versus the user traveled distance. Only signal components that are observable for long traveled distance are shown. The color indicates the estimated received power for the signal components in dB.The RMSE of the user position versus its traveled distance for the outdoor measurements is plotted in Figure11, where the RMSE is averaged over 50 particle filter simulations. As for Figure 8, the red curve denotes the RMSE with no association method applied; the blue and green curves show the RMSE if the ML method and DAS, respectively, are incorporated for data association.Figure 11
The RMSE of the user position versus the user traveled distance for the outdoor measurements. The red curve shows the RMSE if no associations among transmitters are made, the blue curve if the ML method for associations is applied, and the green curve for using DAS.The user is always in a LoS condition to the physical transmitter vTx, and the corresponding LoS signal component is tracked by the user throughout the track, as becomes evident in Figure10. Likewise, the signal component corresponding to the virtual transmitter vTx2 can be tracked after a traveled distance of approximately 22 m until the end.The almost continuous presence of the signals from these two transmitters is reflected in the user RMSE in Figure11. The RMSE without data association methods applied increases in the beginning but then stays constant in the order of 3-4 m with some fluctuations. This is due to the fact that once the variance on the states of transmitters vTx and vTx2 has decreased far enough, they serve as reliable anchors throughout the track. Hence, they prevent the uncertainty about the user state from increasing further, although we measure only the ToA for each signal component in the outdoor measurement scenario.For the same reason, the data association methods cannot really improve the user positioning performance in the outdoor measurements. From another point of view, a correct data association is inherently made for vTx and vTx2 throughout the track, since once these transmitters have been initialized they stay observable throughout the track. However, if a user was to go through the same scenario a second time with prior information of the transmitter states as estimated during the first run, data association would improve the positioning performance as correspondences among transmitters estimated during the first and second run could be found and exploited.In contrast, the user RMSE in Figure8 without data association in the simulations keeps increasing, since there are constantly new transmitters showing up and current transmitters disappear. As mentioned above, the uncertainty about transmitters is high upon initialization, since the measurements obtained from KEST are of fewer dimensions than the transmitter states. In addition, the current user uncertainty adds up to the transmitter uncertainty. No transmitter can be tracked throughout the scenario, and the overall uncertainty keeps increasing. In the simulations, data association relates new transmitters with previously observable transmitters and decreases the uncertainty about the states of transmitters drastically. Consequently, also the uncertainty about and hence the RMSE of the user state decrease.


## 6.1. Simulations in an Urban Scenario

A top view of the urban simulation scenario is depicted in Figure6. The thick black lines represent walls, for example, from buildings, that reflect RF signals, and the black circles are objects such as traffic light poles acting as scatterers. There is one physical transmitter in the scenario marked by the red upward triangle labeled Tx. The user travels with a constant speed of 10 m/s along the blue line with a loop around the central building. The initial and final user positions are labeled START and END, respectively. The traveled distances of the user are marked for every 50 m.Figure 6
The simulation scenario with thick black lines representing reflecting walls and black circles representing objects that scatter the RF signals. The physical transmitter is marked by the red upward triangle labeled Tx. The user travels along the blue line from START to END with one loop around the central building.The transmitter continuously broadcasts a signal that is known to the user and has a rectangular shape in frequency domain with a center frequency of 1.5 GHz and a bandwidth of 100 MHz. As we know the environment, a CIR and the received signal can be modeled for every user position with a simple ray-tracing approach. We incorporate first- and second-order reflections and scattering, that is, single and double reflections and/or scattering. The power loss for the signal being reflected is 3 dB and 6 dB when the signal is scattered at a point scatterer. The average signal-to-noise ratio (SNR) at the user is 7 dB.The user is equipped with an RF receiver and a two-dimensional, rectangular antenna array consisting of nine elements. Hence, both the ToA and the AoA estimates from KEST are incorporated in the estimation of the user and the transmitters’ states. Based on the received signal, KEST estimates the ToAs and AoAs every 50 ms.The results of the KEST estimator are plotted in Figure7. It shows the propagation distance, which is the ToA multiplied by the speed of light, of the signal components versus the traveled distance of the user. Each continuous line represents one signal component and its evolution as the user travels through the scenario. The color of each line indicates the normalized absolute value of the amplitude of the corresponding signal component in linear domain. Since signal components that are observable for a long time can contribute much better to Channel-SLAM than components which are observable only for a short time, only signal components that are observable for a user traveled distance of at least 35 m are plotted and used. Using all detected signal components would dramatically increase the computational complexity and hardly increase the positioning performance.Figure 7
The results of the KEST estimator for the simulations showing the propagation distances of signal components versus the user traveled distance. The propagation distances are the ToA multiplied by the speed of light. Only signal components that are observable for a traveled distance of at least 35 m are shown. The color indicates the normalized amplitude in linear domain.In Channel-SLAM, the user position is estimated relative to the physical and virtual transmitters in the scenario. Thus, to create a local coordinate system, the initial state of the user is assumed to be known. However, no prior knowledge on any transmitter is assumed. In the Rao-Blackwellized particle filter, the number of user particles in the particle filter is 4000, while the number of Gaussian components for each transmitter depends on the first ToA measurement for that transmitter.The root mean square error (RMSE) of the user position versus its traveled distance is plotted in Figure8. The red curve shows the RMSE if no associations among transmitters are made; that is, every signal component that is detected by the KEST algorithm is assumed to be a new transmitter. The RMSEs with the ML association method and DAS from Section 5 being applied are plotted in blue and green, respectively. Since the particle filter is a MC based method, all RMSE curves are averaged over 100 simulations.Figure 8
The RMSE of the user position versus the user traveled distance for the simulations. The red curve shows the RMSE if no associations among transmitters are made, the blue curve if the ML method for associations is applied, and the green curve for using DAS.As we assume the starting position of the user to be known, all three curves start with a low RMSE that increases linearly during the first 200 m as expected. The increase of the RMSE is less due to a bias in the position estimate but more due to an increasing uncertainty, that is, variance, about the user position. After approximately 200 m, the RMSE tends to decrease for all three curves. As more and more transmitters are observed, the weight for some user particles becomes small, and these particles are unlikely to be resampled. Towards the end of the track, the geometrical delusion of precision (GDOP) causes an increase in the RMSE, since most of the transmitters are observed from the same direction. After a traveled distance of around 370 m, several transmitters that had been observed in the beginning are observed again, and correspondences among them can be found. If data association methods are used, the RMSE decreases particularly in that region. The ML method and DAS show a similar performance. Note that there are several reasons for which associations among transmitters can be found. Examples are signal blocking or the geometry of the environment causing virtual transmitters to be observable only from certain regions. In addition, when KEST loses and regains track of a signal component if its received power fluctuates or if another signal component arrives at the receiver with a very small difference in delay, transmitters may be discarded and initialized again at a later point, and associations among them may be found. This explains the increasing positioning performance gain after approx. 50 m using data association.


## 6.2. Outdoor Measurements

In addition to the simulations as described above, we performed outdoor measurements on an airfield. A top view of the measurement scenario is depicted in Figure9. The grey area is an airplane hangar with solid metallic doors. The user track with a total length of 112.5 m is plotted in blue. The user walked along the track starting from the light blue cross labeled START to the black cross labeled END. The traveled distance of the user is marked after 25 m, 75 m, and 100 m. There is one physical transmitter marked by the red upward triangle labeled Tx. The user is in LoS to the physical transmitter throughout the entire track.Figure 9
Top view of the measurement scenario in front of a hangar. The physical transmitter location is marked by the red triangle labeled Tx. The user travels along the blue track from START to END. The expected virtual transmitter locations are marked by the magenta downward triangles.In the scenario, we have three fences labeled Fence 1, Fence 2, and Fence 3. We expect these fences and the hangar door to reflect the RF signal emitted by the physical transmitter. Hence, we expect a virtual transmitter for each of the fences and for the hangar door following Section2.1. The virtual transmitter corresponding to the reflection of the signal at Fence 1 is the magenta downward triangle labeled vTx2. It is located at the physical transmitter position mirrored at Fence 1. Likewise, the location of the virtual transmitter corresponding to Fence 2 is labeled as vTx3. For the reflection of the signal at the hangar doors, we expect the virtual transmitter located at the magenta triangle labeled as vTx1. The expected virtual transmitter corresponding to Fence 3, vTx4, is outside of the boundaries of Figure 9.The Medav RUSK broadband channel sounder [49] was used to perform the measurements. The transmit signal is a multitone signal with a center frequency of 1.51 GHz and a bandwidth of 100 MHz. The signal has 1281 subcarriers with equal gains and a total transmit power of 10 mW.The user was equipped with an RF receiver, recording a snapshot of the received signal every 1.024 ms. For later evaluation, the user carried a prism mounted next to the receiver antenna that was tracked by a tachymeter (Leica Geosystems TCRP1200) to obtain the ground truth of the user location in centimeter accuracy. In addition, the user carried an Xsens MTI-G-700 IMU. Only heading change rate measurements were used from the IMU.On both transmitter and receiver side, single antennas were used. Hence, no AoA information about the impinging signal components can be used for Channel-SLAM. Instead, only ToA estimates from KEST are incorporated. The likelihood function in (19) is adapted accordingly for the evaluation.The results of the KEST estimator for the outdoor measurements are plotted in Figure10. The colors indicate the power estimated by KEST in dBm. As for the simulations, only signal components that are observable for a long user traveled distance are plotted and used. In addition, the ground truth geometrical line-of-sight (GLoS) propagation distances from the physical and the expected virtual transmitters as in Figure 9 to the user are plotted by black lines. They match the KEST estimates very well, justifying the signal model in (5) without considering DMCs in KEST for the measurement scenario.Figure 10
The results of the KEST estimator for the outdoor measurements showing the propagation distances, that is, the ToAs multiplied by the speed of light, of signal components, versus the user traveled distance. Only signal components that are observable for long traveled distance are shown. The color indicates the estimated received power for the signal components in dB.The RMSE of the user position versus its traveled distance for the outdoor measurements is plotted in Figure11, where the RMSE is averaged over 50 particle filter simulations. As for Figure 8, the red curve denotes the RMSE with no association method applied; the blue and green curves show the RMSE if the ML method and DAS, respectively, are incorporated for data association.Figure 11
The RMSE of the user position versus the user traveled distance for the outdoor measurements. The red curve shows the RMSE if no associations among transmitters are made, the blue curve if the ML method for associations is applied, and the green curve for using DAS.The user is always in a LoS condition to the physical transmitter vTx, and the corresponding LoS signal component is tracked by the user throughout the track, as becomes evident in Figure10. Likewise, the signal component corresponding to the virtual transmitter vTx2 can be tracked after a traveled distance of approximately 22 m until the end.The almost continuous presence of the signals from these two transmitters is reflected in the user RMSE in Figure11. The RMSE without data association methods applied increases in the beginning but then stays constant in the order of 3-4 m with some fluctuations. This is due to the fact that once the variance on the states of transmitters vTx and vTx2 has decreased far enough, they serve as reliable anchors throughout the track. Hence, they prevent the uncertainty about the user state from increasing further, although we measure only the ToA for each signal component in the outdoor measurement scenario.For the same reason, the data association methods cannot really improve the user positioning performance in the outdoor measurements. From another point of view, a correct data association is inherently made for vTx and vTx2 throughout the track, since once these transmitters have been initialized they stay observable throughout the track. However, if a user was to go through the same scenario a second time with prior information of the transmitter states as estimated during the first run, data association would improve the positioning performance as correspondences among transmitters estimated during the first and second run could be found and exploited.In contrast, the user RMSE in Figure8 without data association in the simulations keeps increasing, since there are constantly new transmitters showing up and current transmitters disappear. As mentioned above, the uncertainty about transmitters is high upon initialization, since the measurements obtained from KEST are of fewer dimensions than the transmitter states. In addition, the current user uncertainty adds up to the transmitter uncertainty. No transmitter can be tracked throughout the scenario, and the overall uncertainty keeps increasing. In the simulations, data association relates new transmitters with previously observable transmitters and decreases the uncertainty about the states of transmitters drastically. Consequently, also the uncertainty about and hence the RMSE of the user state decrease.


## 7. Conclusion

Within this paper, we derived a novel filtering approach for Channel-SLAM. Using Rao-Blackwellization, the user state is represented by a number of particles and estimated by a particle filter. The states of the landmarks, which are the physical and virtual transmitters in Channel-SLAM, are represented by a sum of Gaussian PDFs, where each Gaussian component PDF is filtered by a UKF. The approach can be applied to SLAM problems in general.We evaluated our approach in simulations in an urban scenario as well as with outdoor measurement data, where we could track a user’s position with only one physical transmitter whose location was unknown. For the simulations in the urban scenario, the user RMSE was always below 21 m. With the presented data association methods applied, it was always below 16.5 m. For the measurements on an airfield, the user RMSE was in the order of 3-4 m.



---

*Source: 4761601-2018-04-23.xml*