Rao-Blackwellized Gaussian Sum Particle Filtering for Multipath Assisted Positioning

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 ourmultipath 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 longterm 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.


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.
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.Section 2 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,  stands for a user particle,  denotes a transmitter or a signal component, ℓ is a component in a Gaussian mixture model, and  stands for a sigma point.
(ii) (⋅)  denotes the transpose of a matrix or vector.
(iii) 1  denotes the identity matrix of dimension  × .
(v) N(x; , C) denotes the PDF of a normal distribution in x with mean  and covariance C.
(vi)  0 denotes the speed of light.

Multipath Assisted Positioning
2.1.Virtual Transmitters.The idea of virtual transmitters is illustrated in Figure 1.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.
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 Reflecting Surface User vTx2 vTx1 vTx3 Scatterer Tx 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).
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 Figure 1, 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.

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 The index  denotes the time instant, the function f  (⋅) is assumed to be known, and k −1 denotes a sample of the process noise with covariance matrix Q.The state is related to the measurement z  by where h  (⋅) is again a known function and n  is a sample of the measurement noise with covariance matrix R. Recursive

Position Estimation Parameter Estimation
Received Signal

IMU Position Estimate
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.
Bayesian estimation works in two steps, the prediction and the update step.The corresponding PDFs can be calculated recursively by for the prediction step and by for the update step, where   is a constant and z 1: denotes the measurements from time instant 1 to .The state transition prior p(x  | x −1 ) and the measurement likelihood p(z  | x  ) are obtained from the movement model in ( 1) and the measurement model in (2), respectively.

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.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 signal (), where the th signal component is defined by a complex amplitude   (  ) and a delay   (  ) at time   .The signal received by the user at time instant   is where () 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 .
The physical transmitter continuously broadcasts the signal () 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 instant  are comprised in the vector where are the ToA estimates for the  TX signal components, or transmitters, and are 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  in  TX .
In the second stage of Channel-SLAM, the user state x u, is estimated simultaneously with the state of the transmitters x TX, .The entire state vector is hence where x ⟨⟩ TX, is the state of the th transmitter.As we consider a two-dimensional scenario, the user state at time instant  is defined by where the user position is defined by p u, = [    ]  and the user velocity by k u, = [V , V , ]  .Each transmitter is defined by its location p TX, = [ TX,  TX, ]  and a clock offset  0, at time instant .The state vector of the th transmitter is hence defined by Our goal is to find the minimum mean square error (MMSE) estimator for x  , which is defined as where z 1: denotes all measurements up to time instant .
We use a recursive Bayesian estimation scheme as in Section 2.2 to estimate the posterior PDF p(x  | z 1: ).This posterior can be factorized as 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 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 the th transmitter is calculated as 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.The heading change rate from the IMU is Δ , .The ToA and AoA measurements for the signal from the th transmitter are  , and  , , respectively, where  , describes the angle between the heading direction k u, of the user and the arriving signal.
With the gyroscope heading change rate Δ , , we predict the movement of the user by where   denotes the time between instants  − 1 and .The two-dimensional rotation matrix R(Δ , ) is defined as where   is the heading noise which is distributed following a von Mises distribution.Hence, the function f  in (1) can be expressed in our case as where the process noise covariance matrix Q is diagonal.As depicted in Figure 3, an AoA measurement for a transmitter  describes the angle  , between the user heading direction k u, 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  2 , and  2 , , respectively, for the th transmitter.Also, we assume no cross-correlation between the single ToA and AoA measurements.The likelihood for the measurement vector z  conditioned on the state vector x  is therefore the product where the predicted ToA between the user and the th transmitter is and the predicted AoA is calculated as The function atan2(, ) calculates the four-quadrant inverse tangent function.It returns the counterclockwise angle between the positive -axis and the point given by the coordinates (, ).

Nonlinear Kalman Filtering
3.1.Unscented Transform.If a random variable x 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 function g(x) of the integration variable x and a Gaussian PDF N(x;  x , C x ) an integration rule of the form can be applied, where X  is the th of the  sig sigma points with its associated weight   .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.

Choice of Sigma Points.
In the literature, different sets of sigma points have been proposed for the unscented transform [44].Let X  be the th sigma point and   its associated weight.The dimension, mean, and covariance of the random variable x are denoted by ,  x , and C x , respectively.In [23], the sigma points and their weights are defined for some  ∈ R as where  = 1, . . ., , (A)  denotes the th row or column of the matrix A, and ( + )C x is factorized into This definition leads to  sig = 2 + 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 .The resulting sigma points are the points in (23) for  = 0. Since the weight of the first sigma point is zero, there are only 2 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.

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 x ⟨⟩ TX, 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  UKF Gaussian components in a Gaussian mixture is estimated by a UKF.The structure of the resulting RBGSPF representation is shown in Figure 4.
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 where x ⟨⟩ u, is the th user particle,  ⟨⟩  its associated weight, and   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 the th transmitter of the th user particle is represented as [48] p (x where z where R ⟨⟩  is the measurement noise covariance matrix for the th transmitter.The number  UKF 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  UKF .The predicted measurement of the ℓth Gaussian component for the th transmitter of the th user particle, ẑ⟨,,ℓ⟩ and the predicted AoA measurement θ⟨,,ℓ⟩ where the th user particle is In the prediction step of the user particle filter, new particles are sampled based on the transition prior p(x u, | x u,−1 ).Hence, the particle x ⟨⟩ u, is drawn as where the function f u, (⋅) describes the user movement model, and k u,−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.

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 the th particle at time instant  is given by [26] This expression can be written as 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 Furthermore, using the Gaussian mixture model representation from (26) where R ⟨⟩  is the measurement noise covariance matrix for the th transmitter and x ⟨,,ℓ⟩ TX, denotes the state of the ℓth Gaussian component of the th transmitter of the th user particle.As we assume no correlation among the ToA and AoA measurements, we have Inserting ( 37) into (36) 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) where the sigma points X  and their weights   can be calculated by (23), where and  is the dimension of a transmitter's state; that is,  = 3.
Finally, the weight of the th particle is calculated as It follows directly from (45)  Note that the weights  ⟨⟩  and  ⟨,,ℓ⟩  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  ⟨⟩  do not depend on the weights  ⟨⟩ −1 from the previous time instant [26].

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  ⟨,,ℓ⟩ −1 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 instant  > 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  UKF .

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 Figure 5, 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).
We define the set of transmitters that had been observed previously but are not detected any more, as old transmitters.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 index  in the association variables for notational brevity.The value of the association variable   denotes an association of the new transmitter with the old transmitter   .We denote the marginalized likelihood of the measurement of the new transmitter that is to be initialized at time instant  by    assuming that the new transmitter is associated with the old transmitter   .The set of association decisions up to time instant  − 1 is denoted by  −1 .From [33], we have where x ⟨,  ⟩ TX, denotes the state vector of the   th transmitter for the th user particle.
Assuming a first-order hidden Markov model, the first integrand in (47) can be simplified to Since we use a Gaussian mixture model to represent the single transmitter states, the second integrand in (47) can be rewritten as Inserting ( 48) and ( 49) into (47) yields which is nonlinear in x ⟨,  ,ℓ⟩ TX, , and rewrite the integral in (50) as Approximating the integral using (22) and inserting it into (50) finally yield The sigma points X  and their weights   can again be calculated by (23) with 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  0 .The set of indices of old transmitters that have not yet been and hence might be associated is denoted by Γ  .
For the ML association method, the association of the new transmitter with the old transmitter   is chosen to be nML, = arg max In DAS, an association is sampled based on the likelihoods    for   ∈ Γ  ∪ {0}, where  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 Algorithm 1.
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.

Evaluations
In the following, we evaluate the RBGSPF derived in Section 4 by means of simulations and actual outdoor measurements.For the evaluations, we implemented a squareroot 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].

Simulations in an Urban Scenario.
A top view of the urban simulation scenario is depicted in Figure 6.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.
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 twodimensional, 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 Figure 7.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.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 Figure 8.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.
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.

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 Figure 9.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.
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 Section 2.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 Figure 10.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.The RMSE of the user position versus its traveled distance for the outdoor measurements is plotted in Figure 11, 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.
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 Figure 10.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 Figure 11.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 Figure 8 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.

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.

Figure 3 :
Figure 3: The user moves in the direction k u, at time instant .The heading change rate from the IMU is Δ , .The ToA and AoA measurements for the signal from the th transmitter are  , and  , , respectively, where  , describes the angle between the heading direction k u, of the user and the arriving signal.

Figure 4 :
Figure 4: Structure of the RBGSPF representation: the user state x u is represented by a number of particles.Each particle estimates the transmitters' states on its own.Each of the  TX transmitter states is represented by a sum of  UKF Gaussian distributions N ℓ with associated weights  ⟨ℓ⟩ .

⟨𝑗⟩ 1 :
are the measurements for the th transmitter and x ⟨,,ℓ⟩ TX,| , P ⟨,,ℓ⟩ | , and  ⟨,,ℓ⟩  are the mean, the covariance matrix, and the weight, respectively, of the ℓth Gaussian component of the Gaussian mixture for the th transmitter.Both x ⟨,,ℓ⟩ TX,| and P ⟨,,ℓ⟩ | are obtained from the update step of the corresponding UKF.Similarly, the likelihood for the measurement of the th transmitter of the th user particle is p (z ⟨⟩  | x ⟨⟩ u, , x and the mean of the corresponding ℓth Gaussian component of the th transmitter is x ⟨,,ℓ⟩ TX,| = [ ⟨,,ℓ⟩ TX,

Figure 5 :
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).

Figure 6 :
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.

Figure 7 :
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.

Figure 8 :
Figure8: 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.

Figure 9 :
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.

Figure 10 :
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.

Figure 11 :
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.