Multipath Assisted Positioning With Transmitter Visibility Information

In multipath assisted positioning, multipath components (MPCs) are regarded as line-of-sight (LoS) signals from virtual transmitters. Instead of trying to mitigate the influence of MPCs, the spatial information contained in MPCs is exploited for localization. The locations of the physical and virtual transmitters are in general unknown but can be estimated with simultaneous localization and mapping (SLAM). Recently, a multipath assisted positioning algorithm named Channel-SLAM for terrestrial radio signals has been introduced. It simultaneously tracks the position of a receiver and maps the locations of physical and virtual radio transmitters. Maps of estimated transmitter locations can be augmented by additional information. Within this paper, we propose to extend the Channel-SLAM algorithm by mapping information about the visibility of transmitters. A physical or virtual transmitter is visible, if its signal is received in a LoS condition. We derive a novel particle filter for Channel-SLAM that estimates and exploits visibility information on transmitters in addition to their locations. We show by means of simulations in an indoor scenario that our novel particle filter improves the positioning performance of Channel-SLAM considerably.


I. INTRODUCTION
Precise localization in indoor and other global navigation satellite system (GNSS) denied scenarios has been both a requirement and an enabler for a variety of services and applications.
For indoor localization, various radio signal based approaches have been introduced, such as fingerprinting [1] or approaches based on the estimation of signal parameters of signal of opportunitys (SoOs) [2], [3]. The SoOs are often network or cellular signals, such as from wireless local area network (WLAN) routers or telecommunication base stations [4], for example. Such signals are available in the majority of scenarios of interest, and no additional infrastructure needs to be deployed to utilize them.
Especially in indoor scenarios, signal propagation delay based positioning methods suffer from multipath and nonline-of-sight (NLoS) propagation [5]. When the transmit signal is reflected and scattered by objects in the environment, The associate editor coordinating the review of this manuscript and approving it for publication was Vittorio Degli-Esposti .
it arrives at the receiver as a superposition of signal components with different delays. A signal component can be the LoS component or a MPC. Due to bandwidth limitations in practical systems, multipath propagation might lead to biased propagation time estimates. Standard techniques try to remove the influence of MPCs on the LoS path [6]. If there is no LoS signal component present, the receiver may untruly deem a MPC to be the LoS signal.
An inherent solution to problems caused by multipath propagation is to exploit the spatial information contained in MPCs. In multipath assisted positioning, each MPC arriving at a receiver is regarded as a LoS signal from a virtual transmitter. Thus, with virtual transmitters, a receiver may be localized with only a single physical transmitter. If the location of the physical transmitter and the geometry of the environment in terms of reflecting and scattering objects is known, the locations of virtual transmitters can be calculated. One of the challenges remaining is the correct association of virtual transmitters to actual MPCs, a problem which is referred to as data association [7], [8]. Various papers have addressed multipath assisted positioning with different technologies and applications such as Long Term Evolution (LTE) [9], ultra-wideband (UWB) [10] or fifth generation (5G) [11], [12] signals, radar [13], or cooperative users [14]. Theoretical bounds were provided in [15], [16].
In a more general case, the environment and the locations of physical transmitters are unknown to the user. Though, the locations of physical and virtual transmitters can be estimated jointly with the receiver position with SLAM [17], [18]. In SLAM terminology, the user equipped with a receiver is localized simultaneously with mapping the transmitter locations as in [19]- [27]. In the following, the term user may refer to the actual receiver depending on the context. The authors of [28], [29] have introduced a multipath assisted positioning algorithm named Channel-SLAM which does not require knowledge about the environment. In a first step, a channel estimator estimates the parameters of signal components and tracks them over time. Such parameters can be the time of arrival (ToA), angle of arrival (AoA), or phase, for example. In a second step, theses estimates are used to obtain the positions of the user and the locations of the transmitters with SLAM. The signal model in Channel-SLAM integrates multiple reflections and/or scattering of the transmit signal at reflecting walls and point scatterers. Geometrical considerations show that the locations of virtual transmitters are static as a user moves through a scenario.
In addition to the locations of transmitters, visibility information, i.e., the information from which positions the LoS signal from a physical or virtual transmitter can be received, can be mapped by a user. In [30], we have shown how visibility information can improve the robustness of data association, i.e, the association among signal components and transmitters, in Channel-SLAM. However, visibility information was not used for estimating the user location. If the user returns to a previously visited location, we expect the same transmitters to be visible as before. Thus, visibility information can be used to improve the positioning performance of Channel-SLAM.
Within this paper, we extend the Channel-SLAM algorithm by mapping and exploiting visibilities of transmitters. A map of transmitter locations created by a user is augmented by the information from where these transmitters are visible. The visibility information can then be used to facilitate the estimation of the user position. We derive a novel particle filter for Channel-SLAM that incorporates visibility information of transmitters by mapping and exploiting this information in a hexagonal grid map. For each hexagon, the probability that a transmitter is visible is estimated. We show with simulations that visibility information does improve the positioning performance of Channel-SLAM considerably.
The remainder of this paper is organized as follows. Section II discusses the idea of multipath assisted positioning and Channel-SLAM. In Section III, we introduce visibility maps and their representation, and derive a novel particle filter for Channel-SLAM exploiting visibility information. Evaluations based on simulations are presented in Section IV. Section V concludes the paper.
Throughout the paper, we use the following notation: • As indices, i identifies a user or user-map particle, j a transmitter or a signal component, and a transmitter particle.
• h denotes the index of a hexagon.
• k is a time instant. • N (x; µ, C) denotes the probability density function (PDF) of a normal distribution in x with mean µ and covariance C.
• c 0 denotes the speed of light. • (·) T denotes the transpose of a matrix or vector.
• · denotes the Euclidean norm of a vector.

II. MULTIPATH ASSISTED POSITIONING A. MULTIPATH ASSISTED POSITIONING AND VIRTUAL TRANSMITTERS
The idea of multipath assisted positioning is depicted in Fig. 1, where one physical transmitter Tx transmits a radio frequency (RF) signal isotropically. For clarity, a possible LoS signal from the physical transmitter to the user is ignored. The transmit signal from the physical transmitter arrives at the user via three different propagation paths. On its first path drawn red, the signal is reflected at the wall, corresponding to a specular reflection of the electromagnetic wave. The user regards the corresponding MPC as a LoS signal from the virtual transmitter vTx1 to the receiver. The location of vTx1 is the location of the physical transmitter mirrored at the wall. The virtual transmitter vTx1 is inherently time synchronized to the physical transmitter.
On the second propagation path drawn brown, the signal is scattered at a point scatterer. A point scatterer distributes the energy of an impinging electromagnetic wave uniformly to all directions. Again, the arriving MPC is interpreted as a LoS signal, now from the virtual transmitter vTx2 to the receiver. The location of this transmitter is at the point scatterer. In the case of scattering, there is an additional propagation distance between the physical transmitter and vTx2, which is the Euclidean distance between the two. Thus, this propagation distance depends on the location of the scatterer. It translates to a time offset τ 0 by dividing by the speed of light c 0 . The time offset can be interpreted as a clock offset or an additional propagation distance between the physical and the virtual transmitter.
On its third propagation path drawn teal, the transmit signal is first scattered at the point scatterer and afterward reflected at the wall. In this case, the location of the resulting virtual transmitter vTx3 is obtained by applying the corresponding two cases iteratively. The location of Tx3 is thus the scatterer location mirrored at the reflecting wall. Accordingly, the time offset between the physical transmitter and vTx3 is again τ 0 .
The third case can be generalized to any number of consecutive reflections and/or scattering. If the transmit signal undergoes only specular reflections at plane surfaces, the corresponding virtual transmitter is time synchronized with the physical transmitter, and accordingly τ 0 = 0 s. If the transmit signal is scattered at least one time, the traveled distance of the signal between the physical transmitter and the last involved scatterer is the additional propagation distance and determines the time offset.
If the signal is only reflected at straight surfaces and scattered at point scatterers, the virtual transmitters are static and independent from the user position. In the sense of virtual transmitters, the effect of diffraction of an electromagnetic wave can be regarded as scattering.

B. SIGNAL MODEL
Within the scope of this paper, we consider a linear and time-variant multipath channel. We assume only one physical transmitter, which is time synchronized to the receiver. The generalization to multiple physical transmitters is straightforward, if the transmit signals are separated in time, frequency or code domain. The transmit signal s(t) from a physical transmitter arrives at the user via different propagation paths and can be modeled at the receiver as a superposition of signal components. At time instant k, the j th signal component is described by a complex amplitude a j (k) and a delay τ j (k), and arrives at the receiver with an AoA θ j (k). Both additive white Gaussian noise (AWGN) and dense multipath components (DMCs) [31] are contained in the colored noise sequence n k (τ ). Considering an antenna array at the receiver, the received signal y m (τ, k) at time instant k at the m th element of the array in the baseband is hence where a m θ j (k) is the response of the m th antenna element. The user samples snapshots of the received signal at time instants k. During the short time interval when the user samples a snapshot, the channel is assumed to be static. While the number of signal components is infinite in theory, it is limited by the sensitivity of the receiver in practical systems.

C. CHANNEL-SLAM
Channel-SLAM exploits the spatial information contained in MPCs for positioning. An overview of the algorithm is depicted in Fig. 2. In the first step, samples of the received signal are used by a channel estimator, which estimates and tracks the parameters of signal components. The results from the channel estimator serve as measurement inputs for the second step, where the user state and the transmitters' states are jointly estimated with SLAM. In the second step, we may use additional sensors such as from an inertial measurement unit (IMU), for example. In addition, prior maps of transmitter locations may be incorporated [32]. In the first step in Channel-SLAM, a channel estimator estimates the parameters of signal components arriving at the receiver. These parameters include the ToA, amplitude or phase, for example. Within the scope of this paper, we use the Kalman Enhanced Super Resolution tracking (KEST) estimator [33]. KEST works in two stages. In the inner stage, the parameters of signal components are estimated and tracked with a Kalman filter. In the update stage of the Kalman filter, a maximum likelihood (ML) estimator such as the Space-Alternating Generalized Expectation-Maximization (SAGE) algorithm [34] is used to obtain estimates of these parameters. The SAGE algorithm is a variant of the Expectation-Maximization (EM) algorithm [35] and estimates the signal parameters jointly. In the outer stage of KEST, the number of signal components is tracked over time by a number of Kalman filters running in parallel. Each Kalman filter carries a different hypothesis on the number of signal components in the received signal. Thus, the estimate of the number of signal components is an additional output of KEST. The tracking nature of the Kalman filters inherently yields associations among signal components for adjacent time instants. Though, KEST cannot find correspondences among signal components that are not tracked continuously. For example, if track of a signal component is lost and the signal component is detected again at a later time instant, KEST cannot associate this newly detected with the original signal component. We denote the number of signal components that KEST tracks at time instant k by N TX,k . The overall number of signal components that have been detected from time instants zero to k and tracked by KEST is denoted by N a TX,k . Every time a new signal component is detected, N a TX,k is increased by one.
There is no differentiation between physical and virtual transmitters in Channel-SLAM, and therefore no differentiation between the LoS component and MPCs of a received signal. Accordingly, the term transmitter comprises both physical and virtual transmitters. Since each signal component corresponds to one transmitter, N TX,k and N a TX,k also denote the number of visible transmitters at time instant k, and transmitters that have been visible during time instants zero to k, respectively.
Within the scope of this paper, we assume that the receivers are equipped with rectangular antenna arrays that are aligned with the user orientation, allowing to exploit the AoA information of signal components. Thus, the ToA estimates τ k and AoA estimates θ k of the signal components estimated by KEST at time instant k are stacked in the radio measurement vector where contains the ToA estimates of the N TX,k detected signal components, and the corresponding AoA estimates. The ToA estimate τ j,k describes the propagation time of the signal traveling from the j th transmitter to the user. The AoA estimate θ j,k is the angle between the user heading and the incoming signal from the j th transmitter.
In the second step of Channel-SLAM, the locations and clock offsets of the transmitters and the user position and velocity are estimated jointly with SLAM. Each signal component detected by KEST corresponds to one transmitter. As described in Section II-A, the locations of the physical and virtual transmitters are static when the transmit signal is reflected at straight walls and/or scattered at point scatterers. Thus, the transmitter model in Channel-SLAM is the same for both the physical and virtual transmitters: each transmitter can be described by a location p TX in two dimensions and an additional propagation distance corresponding to a clock offset τ 0 .
The state vector x <j> TX,k of the j th transmitter at time instant k includes its two-dimensional location p <j> TX,k in Cartesian coordinates and its clock offset τ <j> 0,k , The user state vector x u,k at time instant k consists of the position p u,k = [x k y k ] T and velocity v u,k = v x,k v y,k T of the user in two dimensions each. It is expressed as The entire state vector consists of the user state and the state of the N TX,k transmitters, given by It is clear from the dimension of the user state vector in (6), the transmitters' state vector in (5), and the radio measurement vector in (2) that Channel-SLAM considers an underdetermined system in each snapshot, as there are more unknowns than measurements. Thus, it is crucial that the channel estimator finds correspondences among signal components in consecutive snapshots of the received signal. Only when such correspondences are found, the entire system becomes overdetermined over time due to the static nature of the transmitters. As a consequence, signal components that can be tracked for a long traveled distance of the user can contribute much to Channel-SLAM.
Channel-SLAM seeks to find the minimum mean square error (MMSE) estimatorx 0:k for the user and the transmitter states from time instants 0 to k with SLAM. It is expressed in terms of the posterior PDF p x 0:k |z R,1:k , u 1:k , where z R,1:k denotes the measurements from time instants 1 to k, aŝ The variables u 1:k denote the user control input from time instants 1 to k and will be examined later in this section. The posterior PDF of the state vector can be factorized as p x 0:k |z R,1:k , u 1:k = p x TX,0:k , x u,0:k |z R,1:k , u 1:k = p x u,0:k |z R,1:k , u 1:k × p x TX,0:k |x u,0:k , z R,1:k , (9) separating the user space from the transmitter space. The second term in the last line of (9) denotes the transmitters' posterior PDF conditioned on the user state. We assume no correlation among the estimates of the parameters of signal components, and therefore among the measurements for the single transmitters. The latter assumption does not hold if the parameters of two or more signal components are close to each other in all estimated domains. For example, if the respective differences in the ToA, phase, amplitude and AoA of two signal components are very small, the estimates are correlated. However, if such a situation occurs, we expect it to last for only a very short period of time. Thus, we expect only short term biases in the channel estimator, and unbiased estimates in the long term.
It can be shown [36] that the transmitters' conditional posterior PDF can be written as p x TX,0:k |x u,0:k , z R,1:k = N TX,k j=1 p x <j> TX,0:k |x u,0:k , z R,1:k .
(10) VOLUME 8, 2020 To calculate the posterior PDF, Channel-SLAM applies Bayesian recursive estimation [37], which works in two stages. In the prediction stage, the state evolves following a movement model. In the update stage, the state is updated by measurements, which are the estimates from KEST as in (2).
The transmitters are static in the model of Channel-SLAM. The probabilistic model of the transmitter state evolution is expressed by where v TX,k is a noise sample from a zero-mean distribution with very small variance. While the transmitters are considered static, v TX,k is added for numerical stability.
For the prediction of the user state, we incorporate heading change rate measurements from a gyroscope with the control input u k . The model of the user state evolution is expressed by where the heading change rate measurements are incorporated in the matrix F u,k as in [29], and v u,k is a noise sample from a distribution with a covariance matrix as described in more detail in [38]. It depends on the assumed dynamics of the user, who is assumed to be a pedestrian in Channel-SLAM.
In the update stage, we assume two types of radio measurements, which are ToA and AoA. As mentioned above, they are obtained from the channel estimator in the first stage of Channel-SLAM as in (2).
The likelihood of the measurements from the j th transmitter is calculated by where the predicted measurements are given bŷ for the ToA, and bŷ for the AoA, respectively. The four quadrant inverse tangent function atan2 (y, x) returns the unique counter-clockwise angle between the positive x-axis and the line connecting the origin with the point given by the coordinates (x, y).
The measurement noise samples are assumed to be drawn from zero-mean Gaussian distributed noise sequences. The variance for the ToA measurement is denoted by σ 2 τ,j , and for the AoA measurement by σ 2 θ,j . They are obtained from the covariance matrices of the Kalman filters in the KEST estimator as in [33].
The actual implementation of Channel-SLAM uses a Rao-Blackwellized particle filter to estimate the user and transmitter states jointly. Following (9), a user particle filter tracks the user position and velocity over time. The posterior PDF of the user state is modeled by p x u,0:k |z R,1:k , u 1: where x <i> u,0:k is the i th of the N p user particles, w <i> 0:k its associated weight, and δ (·) the Dirac delta distribution. For each user particle, the transmitter states are estimated independently from other user particles by particle filters as well. The state of the j th transmitter for the i th user particle is therefore expressed by the model where x <i,j, > TX,0:k is the th of the N p,Tx transmitter particles and w <i,j, > 0:k its associated weight. A detailed derivation of the Channel-SLAM algorithm can be found in [29].

III. CHANNEL-SLAM WITH VISIBILITY INFORMATION A. VISIBILITY MAPS
In Channel-SLAM, we map the locations of physical and virtual transmitters. Such maps can be augmented by additional information. As a user moves through a scenario, different physical and/or virtual transmitters are visible from different user locations. We say that a transmitter is visible from a certain user position, if the signal from that transmitter is received by the user in a LoS condition. Following the idea of multipath assisted positioning in Section II-A, a MPC is regarded as the LoS component from a different virtual transmitter. In the rationale of multipath assisted positioning, there is by definition no NLoS propagation. Considering static scenarios, we expect that if a user returns to a previously visited location, the same transmitters are visible as before.
For mapping visibility information, we use a location based map. In particular, visibility information is stored in a hexagonal grid map. The two-dimensional space is discretized with a grid of adjacent hexagons. Each hexagon is assigned a unique index. The entire visibility map M can thus be described as a set of hexagon visibility states, where N H is the number of hexagons. Each hexagon state M h in turn contains information on the visibility of each of the N a TX,k transmitters, The value M <j> h represents probabilistic information on whether the j th transmitter is visible from the h th hexagon or not. We use probabilities since due to the discretization of the space, a transmitter might be visible at one position in the hexagon, but not at another position in the same hexagon. Thus, we have where V <j> h donates the probability that the j th transmitter is visible from an arbitrary position in the h th hexagon, and V <j> h Fig. 3 depicts a simple scenario with a physical transmitter Tx and dark gray lines representing walls that block the signal from the transmitter Tx. Each hexagon is filled according to the probability that the transmitter is visible from an arbitrary user position within that hexagon. In hexagons filled dark, the probability that the transmitter is visible is small, whereas in bright hexagons, the probability is high. In hexagons filled gray with a wall cutting them into two equal halves, the probabilities of the transmitter being visible and being not visible are equal. h,k |x u,0:k , z V,1:k can be estimated, though, where z V,1:k are observations on the transmitters visibilities from time instants one to k and will be covered in more detail later. We assume that each of these PDFs follows a Beta distribution [39]. They represent the belief that a transmitter is visible in a hexagon or not. In general, a Beta distribution with parameters p and q is defined by the PDF for x ∈ (0, 1), and B (x; p, q) = 0 otherwise. The Beta function normalizes the Beta distribution, and it is defined by where (·) is the Gamma function [39]. Fig. 4 shows the PDF of the Beta distribution for different parameters. Intuitively, the Beta distribution as defined above may represent the belief in the value of the random variable x ∈ (0, 1). If the value of the parameter p is high and q is low, the value of x is likely to be high and vice versa. In addition, if the sum of the parameters p and q is high, the degree of trust that we have in our belief of x is high as well. If the sum of the parameters is low, the degree of trust in our belief of x is low. Thus, if the parameters p and q are obtained from observations following a binomial distribution, the number of observations represents how reliable the belief in the random variable x is.
The PDF at time instant k of the j th transmitter being visible in hexagon h is given by the Beta distribution In Section II, we used the ToA and AoA estimates from the inner stage of KEST as measurement inputs in the recursive Bayesian estimation scheme. The correspondences among signal components at adjacent time instants were assumed to be known implicitly. In the following, we drop this assumption, and explicitly introduce an additional measurement vector z V for mapping the visibilities of transmitters. This vector z V is an output of the outer stage of KEST and has N a TX,k entries, where the j th entry at time instant k is denoted by z <j> V,k . If the j th of the N a TX,k transmitters is visible at time instant k, the value of z <j> V,k is the index of the ToA measurement for that transmitter in z R,k . Otherwise, z <j> V,k is zero. Thus, z V,k has N TX,k non-zero entries. Fig. 5 illustrates the relationship of the vectors z R,k and z V,k with an example, assuming for clarity that only ToA measurements are available. The vector z R,k holds the ToA measurements for N TX,k = 3 currently visible transmitters. The transmitters visible at time instant k have the indices 3, 4 and 7, corresponding to the indices of the entries in z V,k that are non-zero. Each of these non-zero entries denotes the index of the ToA measurement for the respective transmitter in z R,k . For example, z <3> V,k = 1, meaning that the ToA measurement of the third transmitter is the first entry in z R,k . In total, N a TX,k = 8 transmitters have been visible, i.e., N a TX,k = 8 signal components have been tracked by KEST up to time instant k. On the one hand, z V,k can be regarded as an association vector, as it associates signal components at time instant k with signal components from previous time steps. On the other hand, it can be regarded as a measurement for visibilities of transmitters, as it indicates which transmitters are visible at which time instant. Hence, the vector z V,k may be named visibility measurement vector in the following. The overall measurement vector z k comprises the estimates from both the inner and outer stage of KEST. Thus, we define where z R,k is defined as in (2). It contains the ToA and AoA estimates from the inner stage of KEST for the currently visible N TX,k transmitters.

C. BAYESIAN FORMULATION AND RAO-BLACKWELLIZED PARTICLE FILTER FOR CHANNEL-SLAM
We seek to estimate the states of the user, x u,0:k , the locations of the transmitters x TX,0:k , and the visibility map M 0:k given the measurements z 1:k and the user control input u 1:k with recursive Bayesian estimation. The dynamic Bayesian network (DBN) in Fig. 6 shows the dependencies of the involved variables. We assume a first-order hidden Markov model. In comparison to the standard Channel-SLAM algorithm, the visibility measurements z V and map M are added. The corresponding posterior PDF can be written as p (x 0:k , M 0:k |z 1:k , u 1:k ) = p x TX,0:k , x u,0:k , M 0:k |z 1:k , u 1:k = p x u,0:k , M 0:k |z 1:k , u 1:k × p x TX,0:k |x u,0:k , M 0:k , z 1:k , u 1:k = p x u,0:k , M 0:k |z 1:k , u 1:k p x TX,0:k |x u,0:k , z 1:k , (27) where we assume that the transmitter states are conditionally independent from the user control input u 1:k and the visibility map M 0:k . The latter assumption is based on that without knowledge on the geometry of the environment, a virtual transmitter's location does not reveal much information about its visibility. For example, in Fig. 1, the transmitters Tx and vTx2 are visible in the vicinity of the corresponding transmitter's location, while vTx1 and vTx3 are not. In addition, we do not differentiate between physical and virtual transmitters in our model, and assume no knowledge on the propagation path of a signal component.
Given the structure of the last line of (27), we use a Rao-Blackwellized particle filter [40] to jointly estimate the user state and the visibility map, and to simultaneously estimate the transmitter locations. We name the particle filter estimating the user state and the visibility map the user-map particle filter, and the particle filters estimating the transmitter states transmitter particle filters. The history of the i th user-map particle at time instant k is of the form and has an associated weight w <i> 0:k . The posterior PDF in the user-map particle filter, i.e., the first factor in the last line of (27)

D. DERIVATION OF THE WEIGHTS FOR THE USER-MAP PARTICLE FILTER
In general, it is hard to draw samples from the PDF on the left-hand side of (29). Instead, we use the idea of importance sampling. This idea is to draw samples from an importance density q x u,0:k , M 0:k |z 1:k , u 1:k , and compensate for the mismatch by adapting the weights on the right-hand side of (29). If we can evaluate the left-hand side of (29), the weights w <i> The general derivation of (32) can be found in [41], and a derivation for the user-map particle filter is presented in Appendix I. With (26), the second factor in the second to last line in (32) can be factorized as p z k |x u,0:k , M 0:k , z 1:k−1 , u 1:k = p z V,k |x u,0:k , M 0:k , z 1:k−1 , u 1:k × p z R,k |x u,0:k , M 0:k , z V,1:k , z R,1:k−1 , u 1:k = p z V,k |x u,k , M k p z R,k |x u,0:k , z V,1:k , z R,1:k−1 .
In the last line of (33), the measurement z V,k depends only on the current map and the user location. In fact, the user state is needed only to indicate the hexagon in which the user is located at time instant k. The radio measurement z R,k is independent from the map and the control input. For the numerator of the last line of (32), we obtain In general, the user position is conditionally dependent on the estimate of the visibility map. However, without a measurement z V of the visibility, the information in the visibility map cannot be used, since the association of the currently visible transmitters to transmitters in the map is not known. Thus, in the first factor in the last line of (34), the current user state x u,k is independent from the visibility map M k−1 .
In the second factor in the last line of (34), the information of the user state x u,k−1 that is relevant for M k is already contained in M k−1 . Thus, M k is independent from the user state at time instant k − 1. In addition, the visibility map is independent from any control input u k .
A crucial step in a particle filter is the choice of the importance density in the denominator of (32), as new particles are sampled from it. The importance density can be written as Following the structure of (35), we define the importance density such that q x u,k , M k |x u,0:k−1 , M 0:k−1 , z 1:k , u 1: where the first term on the right hand side of (36) is the a-priori PDF for the user. For sampling a new user-map particle, the importance density q x <i> u,k , M <i> k |x <i> u,0:k−1 , M <i> 0:k−1 , z 1:k , u 1:k as in (36) is evaluated from the left to the right. First, a new user particle x <i> u,k is sampled with the user movement model in (12) given the old particle x <i> u,k−1 . Then, a new visibility map is sampled based on the new user particle. To sample a new visibility map, we proceed as follows. As the user travels through the scenario, they store counts on how often a transmitter was visible or not in each hexagon for each particle. At time instant k, the j th transmitter has been visible C <i,j> h,k times in hexagon h for the i th user-map particle, and it has been not visibleC h,k−1 are updated for the current hexagon h, in which the user-map particle is located, for each transmitter depending on whether the transmitter is visible or not.
The Beta distribution is the conjugate prior of the binomial distribution. Updating a prior Beta distribution with observations following a binomial distribution results again in a Beta distribution with new parameters. The observation if a transmitter is visible in a hexagon at one time instant may be regarded as a realization of a random variable following a binomial distribution. Thus, the parameters of the Beta distribution representing the belief of the i th user-map particle that the j th transmitter is visible in the h th hexagon are updated accordingly to The denominator in the last line of (38) can be calculated with Bayes' theorem as Thus, the weights become The visibility measurement z V,k is only dependent on the hexagon the user-map particle is in at time instant k. We identify this hexagon by the index h. In addition, the visibilities for the transmitters are assumed to be independent from each other. The factor in the last line of (40) is therefore where M <i,j> h,k−1 refers to the the visibility of the j th transmitter for the i th user-map particle in the hexagon with index h at time instant k − 1.
The likelihoods with respect to the visibility for the j th transmitter, i.e., the factors in the product in (41), can be calculated as the expectation value of the belief that the transmitter is visible or not, The likelihoods with respect to the radio measurements are p z R,k |x u,0:k , M 0:k , z V,1:k , z R,1:k−1 , where we set the likelihood p z R,k |z V,k , x <i> u,k , x <i,j, > TX,k for the j th transmitter to one if the j th entry in z V,k is zero. If that entry in z V,k is non-zero, the likelihood can be calculated following (13), where the indices for the transmitters need to be adapted. A derivation of (43) following [29] can be found in Appendix II. Finally, plugging (41) and (43) into (40), the weights can be calculated as From (44) follows the weight update for the transmitters. The weight of the th particle in the transmitter particle filter for the j th transmitter and the i th user-map particle is The increase of computational complexity in the particle filter derived above compared to the standard Channel-SLAM particle filter is limited to evaluating the terms in (42). For each user-map particle and each hexagon visited by the usermap particle, the two parameters of a Beta distribution need to be stored for every transmitter. However, information on the visibility of transmitters may be used not only in the particle filter, but also for map matching when maps are exchanged, and for data association of transmitters as in [30].

A. SIMULATION SCENARIO
A top view of a mall serving as the indoor simulation scenario is depicted in Fig. 7. In the scenario, there is one physical transmitter, which is depicted by the red triangle and labeled Tx. It continuously transmits a signal which is known to the receiver and has a bandwidth of 100 MHz. The center frequency of the signal is at 1.9 GHz. The signal energy is spread uniformly across the signal bandwidth around the center frequency. The thick black lines in the scenario represent walls which reflect the transmit signal. Likewise, the thick black dots model point scatterers which spread the energy of an impinging signal uniformly to all directions. Based on the locations of the physical transmitter and the walls and scatterers, we can calculate a channel impulse response (CIR) for every user location with ray-tracing. For the simulations, we incorporate virtual transmitters of orders one and two, corresponding to single and double reflections and/or scattering of the transmit signal. For each user location, the CIR is convolved with the transmit signal and AWGN is added to create the received signal.
The update rate of the receiver is 10 Hz. A sampled snapshot of the received signal is thus fed to the KEST estimator every 100 ms. We also use a gyroscope at the receiver. The gyroscope is incorporated as control input in the user movement model in (12). For the evaluations, we assume the starting location and velocity of the user to be known within their local coordinate system. The locations of the walls, scatterers or any transmitters are unknown to the user.
A large hexagon size in the visibility map leads to a large discretization error. If the hexagon size is very small, it becomes unlikely that a user revisits hexagons, and the memory requirement for the map is large. We have found in our simulations that a hexagon side length of 2 m is a good tradeoff and set the size accordingly.
The number of user-map particles in the simulations is 800. The number of transmitter particles is different for every user or user-map particle, for every time instant and for every transmitter. When a new transmitter is initialized, the number of transmitter particles is based on the radio measurement. This number is subsequently adapted to the uncertainty about the transmitter's state with a grid based particle reduction method from [42] as the user traverses through the scenario. The lower the uncertainty about a transmitter state, the less particles are used. For the first user track depicted in blue in Fig. 7, the maximum number of transmitter particles we have observed for one transmitter and one user-map particle was 864. The mean of transmitter particles for that transmitter while it was visible was 492. The physical transmitter in the first user track is initialized with 89 particles for each user-map particle. Averaging over the user-map particles, this number drops down to 11 after a traveled distance of approximately 77 m, and stays at approximately 10 after a traveled distance of approximately 139 m.
During the first user track, 47 different transmitters were detected by KEST. For the other four tracks, the overall number of transmitters detected by KEST were 34, 23, 33, and 30, respectively.

B. SIMULATION EVALUATIONS
The simulation results are depicted in Fig. 8. For each of the five tracks, the positioning error, namely the mean absolute error (MAE), of the user is plotted by the solid lines versus the traveled distance for two cases. The curves in blue, denoted by 'Without Visibility Mapping', are the MAEs of the standard Channel-SLAM algorithm if no visibility mapping is performed. For the curves in red, denoted by 'With Visibility Mapping', the new particle filter derived in Section III is used. It incorporates the visibility information from a prior VOLUME 8, 2020 transmitter map that is created based on the scenario. We use a prior map to exploit the full potential of visibility mapping. This prior visibility map could stem from a different user and is structured in hexagons as described in Section III-A. If the j th transmitter in the map is visible from the center of the h th hexagon, the parameters of the corresponding Beta distributions in the prior map are set to ν The prior map comprises visibility information for all transmitters up to an order of two. For associating signal components detected by the KEST estimator in Channel-SLAM with transmitters from the prior map, the data association method from [30] is used.
In addition to the MAEs, the 95 th and 5 th percentiles are plotted for each track by the dotted and dashed lines, respectively. Hence, the MAE was below the dotted lines in 95% of the simulation runs, and below the dashed lines in 5% of the runs.
Since the starting location and velocity are assumed to be known, all MAE curves in Fig. 8 start at an error of zero at the first time instant, and increase. The error with visibility mapping is then in almost all tracks considerably below the error without visibility mapping. The increasing error near the end of Track 1, Track 3 and Track 5 can be explained with a high geometrical dilution of precision (GDoP). At these positions, all detected signal components arrive at the user from a similar direction.
All error curves in Fig. 8 are averaged over 250 simulation runs.
The MAEs at the end of the first track are in the order of four to ten meters. For comparison, we also performed simulations where a user does not perform Channel-SLAM, but estimates their state only based on the gyroscope and the movement model. I.e., no radio signals are used in this case. Even if the velocity is known perfectly at the beginning of the track, the corresponding MAE grows quadratically and finally is in the order of 320 m for the first track. The positioning performance is thus dominated by the movement model only for a short distance at the beginning of a run, when the uncertainty for all transmitters is high. Once the transmitter state estimates have converged, the influence of the movement model is small.
The evaluations above consider an approximation of a converged prior map of transmitter visibilities that is known to the user to illustrate solely the effect of mapping of and using information about transmitter visibilities. In reality, such a prior map is created by users with crowdsourcing. A first user maps the transmitter states and visibilities with Channel-SLAM and hands the map over to a second user. Since the second user does in general not know their starting location, they need to estimate the translation and rotation relating the coordinate system of the two users in order to be able to use the map. This estimation is referred to as map matching and can be done based on transmitter states estimated by the second user and in the map as in [43]. Once a map match has been found, the second user exploits the information in the map, updates it with own observations. The map is then handed over to a third user, who again performs map matching, uses and updates the information in the map, hands the map on to the fourth user, and so on. We refer to this approach as collaborative Channel-SLAM.
We have evaluated collaborative Channel-SLAM for ten different users walking on tracks that are plotted in Fig. 9. The start and end positions of the u th track are labeled Su and Eu, respectively. The MAEs for five different users averaged over 250 simulation runs are exemplarily plotted in Fig. 10 for two cases with the solid lines. The u th track is labeled Track cu. In the first case represented by the blue lines, no visibility mapping is performed. In the second case represented by the red lines, the visibilities of transmitters are mapped and exploited as in Section III. The dotted and dashed lines represent the 95 th and 5 th percentiles, respectively.
The positioning errors for the cases where visibility information is used is in almost all cases below the error for the case without visibility. In the collaborative approach, visibility information on the one hand improves the positioning performance directly as in Fig. 8. On the other hand, a better positioning performance leads to a more accurate and robust map matching, which in turn improves the positioning performance as well.
The 5 th percentiles are very similar for most of the tracks in Fig. 10, while the 95 th percentiles are most of the time lower if visibility information is used. We see in the simulations that user particles with a large bias in the position estimate are resampled faster in the particle filter. This leads to the conclusion that visibility information mainly improves the robustness of Channel-SLAM by preventing large biases in the user position estimates.

V. CONCLUSION
Channel-SLAM is a multipath assisted positioning SLAM algorithm. The position of the user is estimated jointly with building a map of physical and virtual transmitter locations based on the estimated parameters of signal components. Within this paper, we have augmented such a transmitter location map by the information from which locations a transmitter is visible. Information on transmitter visibilities is mapped in a hexagonal grid map. For each hexagon, the probability that a transmitter is visible is represented by a Beta distribution. We have derived a new particle filter for Channel-SLAM that incorporates both the signal components' estimated parameters as well as estimates on transmitter visibilities. If a user revisits a hexagon, or if a prior map with transmitter visibilities is available, the visibility information is used to update the weights in the particle filter. For transmitter visibilities, no new measurements as such are required. Instead, the visibilities are estimated based on information that is already available from the channel estimator.
Our simulations show that if a converged map of visibility regions is available as prior information, it increases the positioning performance of Channel-SLAM for a single user considerably. In the case of collaborative Channel-SLAM, where a map of visibility regions is estimated by different users subsequently, visibility information additionally increases the positioning performance by a higher map matching robustness.

APPENDIX I. WEIGHT DERIVATION
In the following, we derive (32) from (30). The numerator in (30) where we make use of the fact that the current user-map particle does not depend on previous measurements and control inputs if its state at time instant k − 1 is known. In addition, we do not incorporate information from future time instants.
Inserting (47) exploiting that the denominator in the second line of (46) is not dependent on the user or the visibility map. The weights in (32) result from inserting (48) and (31) into (30).

APPENDIX II. RADIO MEASUREMENT LIKELIHOOD
In this section, the likelihood in (43) is derived. First, we marginalize over the transmitter state at time instant k. VOLUME 8, 2020 With the assumption of independent transmitter measurements the radio likelihood is p z R,k |x <i> u,0:k , z V,1:k , z R,1:k−1 = N a TX,k j=1 p z R,k |z V,1:k , z R,1:k−1 , x <i> u,0:k , x <i,j> TX,k × p x <i,j> TX,k |x <i> u,0:k , z V,1:k , z R,1:k−1 dx <i,j> TX,k . (49) The first factor in the integral of (49) is p z R,k |z V,1:k , z R,1:k−1 , x <i> u,0:k , x <i,j> TX,k = p z R,k |z V,k , x <i> u,k , x <i,j> TX,k . (50) The radio likelihood z R,k is conditionally dependent on the visibility measurement z V,k for a correct association of the radio measurements to the transmitters. It is independent from any previous states or measurements. The second factor in the integral of (49) The last term in the last line of (52) is obtained from the movement model of the transmitter given by (11). Thus, inserting (50) and (52)