Positioning Using Terrestrial Multipath Signals and Inertial Sensors

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


Introduction
Today, most smartphones are equipped with global navigation satellite systems (GNSSs) receivers which allow using applications on the smartphones for navigation [1].GNSSs provide sufficient position accuracies for mass market application in open sky conditions.However, indoors or in urban canyons the GNSS positioning accuracy could be drastically reduced.In these situations, the GNSS signals might be blocked, degraded by multipath effects, or received with low power.To enhance the positioning performance indoors, different methods and sensor systems can provide position information rather than relying on GNSSs [2][3][4].Most of the indoor positioning systems use local infrastructure like positioning with Radio Frequency Identification (RFID) [5], mobile communication base-stations [6,7], wireless local area network (WLAN) [8], or ultra-wideband (UWB) [9][10][11].However, also these wireless radio technologies experience multipath and non-line-of-sight (NLoS) propagation.
Multipath propagation is experienced when the transmitted signal arrives at the receiver via several propagation paths.These propagation paths with different delays are caused by reflections, diffractions, and scattering of the electromagnetic wave.Hence, the signal at the receiving antenna consists of a superposition of multiple replicas of the transmitted signal, where each version is called multipath component (MPC) traveling along an individual propagation path.The delay estimate of standard algorithms like the delay locked loop (DLL) is biased in multipath propagation environments [12].Algorithms like [13][14][15] reduce the multipath error by modifying the DLL structure.Other algorithms estimate the channel impulse response (CIR) in order to mitigate the influence of multipath propagation on the delay estimate, for example, [16][17][18][19][20]. To retrieve the required delay from the CIR, the path with the smallest delay is treated as the line-of-sight (LoS) path.However, treating the smallest delay as the LoS path may result in weak positioning performance in NLoS situations.Furthermore, even advanced multipath mitigation algorithms reduce the multipath effects only to a certain degree due to limited signal bandwidth and measurement noise [18].
Nowadays, multipath exploitation instead of mitigation is attracting more and more interest.The authors of [21,22] exploit multipath propagation for positioning of mobile terminals using multipath fingerprinting algorithms.Other algorithms, for example, [23,24], interpret reflected signals as signals emitted from virtual transmitters (VTs), where the VT positions are precalculated based on the knowledge of the reflecting surface and physical transmitter positions.Furthermore, the authors of [25] estimate and track the phase information of MPCs using an extended Kalman filter (EKF) and estimate the user position using a time difference of arrival (TDOA) positioning approach.Other algorithms like [26] use a nonlinear least squares algorithm combining UWB measurements at several receiver positions to estimate the positions of the VTs and the receiver simultaneously within small scale scenarios.
This paper describes and extends the multipath assisted positioning algorithm referred to as Channel-SLAM; see [27][28][29][30][31]. Channel-SLAM considers a moving receiver and is suitable for GNSS denied areas like indoor areas.Similarly to other multipath assisted positioning approaches, Channel-SLAM interprets MPCs as LoS signals emitted from VTs.In addition to reflected signals, Channel-SLAM considers also paths occurring due to multiple number of reflections, diffractions, or scattering as well as combinations of these effects.As a consequence, the reception of several MPCs allows position estimation even if only one physical transmitter is present.Interpreting MPCs as directly propagated signals originated from VTs, Channel-SLAM enables positioning also in NLoS situations.Additionally, Channel-SLAM does not require any prior knowledge on locations of reflecting surfaces as Channel-SLAM estimates the receiver position, velocity, clock bias, and the VT positions simultaneously which can be interpreted as simultaneous localization and mapping (SLAM) with radio signals.In [27,28,31], we showed that positioning is possible in NLoS scenarios using MPCs without the knowledge of the room geometry by using Channel-SLAM.We investigated in [27] TDOA positioning and especially TDOA between MPCs such that time synchronization between physical transmitters is not essential.In [31], we derived Channel-SLAM based on a Rao-Blackwellized particle filter (RBPF) and compared the accuracy of Channel-SLAM to a derived posterior Cramér-Rao lower bound (PCRLB).However, the Channel-SLAM algorithms in [27,28,31] use linear antenna arrays and assume the knowledge of the physical transmitter position.
In this paper, we propose an implementation of Channel-SLAM that uses only a single receiving antenna and fuses similarly to [29,30] additional information obtained from an inertial measurement unit (IMU).Today many smartphones feature Microelectromechanical System (MEMS) IMUs, which can provide short term relative orientation and position information.Theoretically, the measurements of the IMU can be directly used in an inertial navigation system.However, the position calculation involves double integrations; hence, even small measurement errors quickly cause a drift in the position solution [32].To avoid that, we only fuse heading measurements from the IMU which solely requires an alignment of the coordinate systems.The heading information of the IMU allows improving the performance of Channel-SLAM by resolving ambiguities and angle of arrival (AoA) measurements are not mandatory anymore.Being a relative positioning system, Channel-SLAM requires an initial prior knowledge of the receiver position and moving direction to define the coordinate system.The positioning algorithm derived in this paper is based on a RBPF where we employ a new transition model for pedestrians.In [29,30], we showed that positioning with only one physical transmitter is possible if MPCs and heading information from an IMU are used.Compared to [29,30], the novel transition model enables a performance gain in the position accuracy.In addition to [27][28][29][30][31], we propose a method to dynamically adapt the number of particles which significantly reduces the computational complexity.Furthermore, a PCRLB for Channel-SLAM is derived which incorporates heading information obtained by using an IMU.The developed positioning algorithm is evaluated based on measurement data obtained in an outdoor scenario, where the position of the physical transmitter is unknown.Based on these measurements, we compare the accuracy of Channel-SLAM to that of the derived PCRLB.
The paper is structured as follows: Section 2 describes the signal model; afterwards, Section 3 describes the proposed algorithm which is split into four subsections: Section 3.1 addresses Channel-SLAM; Section 3.2 describes two different transition models using the heading information from an IMU; Section 3.3 summarizes the RBPF; Section 3.4 describes the implementation of the RBPF; afterwards, we derive in Section 4 the PCRLB for Channel-SLAM incorporating the heading changes of the IMU.Thereafter, Section 5 evaluates the algorithm based on measurement data.The last section, Section 6, concludes the paper.
Throughout the paper, we will use the following notations: (i) [⋅]  stands for the vector transpose.
(ii) All vectors are interpreted as column vectors.
(iii) Vectors are denoted by bold small letters.
(ix) p() denotes the probability density function of .
(x)  is the speed of light.
(xi) x denotes the estimation of .

Concept of Virtual Transmitters
Mathematically, the behavior of the multipath channel can be described by the time variant CIR ℎ(  , ), where   indicates the discrete time instants and  the delay [33].According to [33], the CIR ℎ(  , ) can be assumed to be constant for a short time interval  at discrete time   with index , for  0 ≤   ≤  0 + , where (  ) is the number of MPCs,   (  ) is the delay,   (  ) the complex amplitude of the th MPC, and () stands for the Dirac distribution [34] (please note that the CIR is generally a summation of an infinite number of MPCs; however, a practical receiver is only capable of capturing signals whose powers are above a certain sensitivity level).For notational conveniences, the LoS propagation path is considered also as a MPC in this paper.
Assuming that the transmitted signal (  ) is band-limited with bandwidth  and time-limited with a length smaller than , the signal received at time   sampled with rate , bin indices  = 0, . . .,  − 1, and the delay   = / can be expressed as where (  ) denotes the white circular symmetric normal distributed receiver noise with variance  2  .Using vector notation we obtain from (2) y (  ) = [ (  ,  0 ) , . . .,  (  ,   ) , . . .,  (  ,  −1 )] .(3) In order to obtain the sparse structure of the CIR from the measurements y(  ), super resolution multipath estimation algorithms are necessary.The received signal is geometrically dependent on the transmitter and receiver positions as well as on the environment.Thus, the channel is spatially correlated as long as the spatial sampling is small enough.Hence, we use in this paper the dynamic multipath estimator named Kalman enhanced super resolution tracking (KEST) [20,[35][36][37] for estimating and tracking multipath parameters.KEST allows estimating the evolution of the CIR over time which is essential for Channel-SLAM as shown in the following section.KEST consists of a Kalman filter (KF) to estimate the complex amplitude α (  ) and delay τ (  ) for each MPC  utilizing maximum likelihood (ML) estimates as measurements.In the used implementation, KEST uses a standard model for the CIR which comprises a sum of weighted Dirac impulses as in (1).This model describes distinct paths sufficiently well.However, dense multipath components (DMC) lead to a model mismatch in the used KEST implementation.This model mismatch results in an increased variance of the estimated MPC parameters used as measurement noise in Channel-SLAM.For further details about KEST, see [20,[35][36][37].
To use the delay measurements of the tracked MPCs for positioning, a model describing the delays   (  ) depending on the current user position r  (  ) is necessary.For developing such a model, we consider a static environment with a fixed transmitter and a receiver moving along an arbitrary trajectory.Figure 1 summarizes four propagation scenarios; for a detailed description see [31].In the first scenario, the transmitted signal is reflected on a reflecting surface indicated by the blue lines.For reflection, we consider the effect of an electromagnetic wave reflected by a reflecting surface.When the receiver is moving, the reflection point on the reflecting surface is moving as well.If we mirror the physical transmitter position on the reflecting surface, we obtain the position r VT,1 of VT 1 which is static during the receiver movement.The distance between VT 1 and the receiver is equivalent to the propagation time of the reflected signal multiplied by the speed of light.Hence, the reflected signal can be interpreted as a direct signal from VT 1 to the receiver.
This behavior can be extended to a multiple reflection scenario represented by the red lines.The transmitted signal is reflected two times.Equivalently, the location of VT 2 can be determined by mirroring the transmitter position at both reflecting surfaces, as indicated in Figure 1.The distance between VT 2 and the receiver is equivalent to the propagation time of the reflected signal multiplied by the speed of light.Thus, the signal reflected twice can also be interpreted as a direct signal from VT 2 to the receiver.
Figure 1 exploits by the orange lines additionally a scenario where the signal is scattered, for example, at a lamp post.The propagation effect of scattering occurs if an electromagnetic wave impinges on an object and the energy is spread out in all directions [38].Geometrically, the effect of scattering can be described as a fixed point  at position r  in the pathway of the MPC.We define  as VT 3 at the position r  which is constant for all receiver positions for the MPC.Additionally, we treat  VT > 0, the constant distance between physical transmitter and scatterer, as an additional propagation distance associated with the MPC.Hence, the scattered signal can be interpreted as a direct signal from VT 3 to the receiver, however, with a constant offset  VT .Scattering and diffraction can be geometrically described as a fixed point  at position r  in the pathway of the MPC and are considered as one model.Hence, unless otherwise stated, the description of scattering is equivalent for diffraction.
The fourth scenario considers the combination of both effects indicated in green.The transmitted signal is scattered at  and afterwards reflected on the first reflecting surface.When the receiver is moving, the reflection point on the reflecting surface is moving as well.Hence, VT 4 is defined by mirroring the scatterer  at the first reflecting surface.Furthermore, between the transmitter and  additional interactions are possible leading to the same position of VT 4 .
To summarize, the propagation path of the th MPC can be equivalently described as a direct path with propagation The figure shows four propagation scenarios: First scenario (blue): the transmitted signal is reflected on a reflecting surface.VT 1 is defined by mirroring the physical transmitter position at the surface.Second scenario (red): the transmitted signal is reflected twice.VT 2 is defined by mirroring the physical transmitter position at both surfaces.Third scenario (orange): the transmitted signal is scattered at . VT 3 is defined at the position of .Fourth scenario (green): the transmitted signal is scattered and afterwards reflected on a reflecting surface.VT 4 is defined by mirroring the scatterer  at the surface.
length   (  ) between VT  and the receiver plus an additional constant propagation length  VT, (  ); hence, where  denotes the speed of light and r VT, (  ) the position of the th VT (please note that the position of the VTs and the additional propagation lengths are constant over time.Nevertheless for notational convenience a time dependence on   is introduced here).The additional propagation length is zero, that is,  VT, (  ) = 0, if only reflections occurred on the pathway between physical transmitter and receiver or greater than zero, that is,  VT, (  ) > 0, if the MPC is interacting with at least one scatterer.In general,  VT, (  )/ can be interpreted as a clock offset between the th VT and the physical transmitter.

Channel-SLAM
in Channel-SLAM with the corresponding variances   (  ).
Because the VT positions are unknown, the receiver position and the positions of the VTs have to be estimated simultaneously.Thus, the state vector x(  ) describing the complete system at time instant   for (  ) MPCs is with the receiver states x  (  ) and the VT states x VT (  ).The receiver state x  (  ) includes the receiver position r  (  ), the receiver velocity k  (  ), and the receiver's clock bias   (  ); hence, According to the description given in the previous section and ( 4), an MPC can be represented by a direct path between a VT and the receiver plus an additional propagation length.Hence, the parameters representing the th VT are defined as where r VT, (  ) is the position of the th VT and  VT, (  ) the additional propagation length.Using vector notation for all VTs, we obtain Additionally, as illustrated in Figure 2, an IMU is used.The IMU provides measurements of the acceleration a  (  ) and turn rates    (  ) in three dimensions.After calibration,

Receiver IMU KEST Calibration
Channel-SLAM x u (t k−1 ) x u (t k ) x ６４ (t k−1 ) the heading change Ψ(  ) is used in Channel-SLAM as a control input and is therefore directly integrated into the transition model.
We use a discrete time representation for the transition and measurement model of the dynamic system with The transition model in (10) describes the state evolution from time instant  −1 to time instant   employing a possible nonlinear function f(⋅, ⋅, ⋅) with the process noise n  (  ) and using a control input which is in our case the heading change Ψ(  ).The control input is considered as perfectly known and hence error-free.The measurement model (11) relates the state vector to the measurements by a possible nonlinear function h(⋅, ⋅) and the measurement noise n ℎ (  ) at time instant   .Figure 3 shows the considered dynamic Bayesian network, that is, a first-order hidden Markov model.Equations ( 10) and ( 11) can also be interpreted from a Bayesian perspective: based on measurements, we want to recursively estimate the unknown probability density function (PDF) of the state x(  ).In a recursive Bayesian formulation, this problem can be described as finding the posterior probability distribution Recursive Bayesian filtering provides a methodology to optimally estimate ( 12) by a prediction step to calculate p(x(  ) | z( 1:−1 ), Ψ( 1: ), x( 0 )) and an update step to obtain p(x(  ) | z( 1: ), Ψ( 1: ), x( 0 )) which considers the measurement z(  ) at time instant   with the likelihood function p(z(  ) | x(  )) [39,40].By assuming independence between the transition priors of the receiver state vector x  (  ) and the VT state vectors x VT, (  ) associated with the MPCs  = 0, . . ., (  )−1, the transition prior p(x(  ) | x( −1 ), Ψ(  ), x( 0 )) is defined here as where we inherently assume independence among MPCs, that is, propagation paths interacting with distinct objects.This is based on the well-known uncorrelated scattering assumption in wireless propagation channel modelling [38].
We obtain for the transition prior p(x VT, ( For the transition prior p(x  (  ) | x  ( −1 ), Ψ(  )) of the receiver state vector we provide in Section 3.2 two models indicated by the function f(x  ( −1 ), Ψ(  ), n  (  )) in Figure 3.
Assuming the elements of z(  ) to be independent Gaussian distributed conditioned on the current state x(  ), p(z(  ) | x(  )) can be expressed as with the propagation length for the th MPC, where  2 , (  ) denotes the corresponding variances.

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

Gaussian-Transition-Model.
The first proposed transition model is based on a discrete white noise acceleration model [47], referred to as Gaussian-Transition-Model, with and the velocities where V , (  ), V , (  ) are the corresponding velocities in - direction and the receiver's clock bias   (  ) where a standard clock bias model is used [12,48] (Please note that if transmitter and receiver oscillators provide different frequencies, a clock drift parameter has to be considered additionally).The transition matrix A  (  , Ψ(  )) in ( 17) includes a rotation matrix with the heading changes Ψ(  ), with where   =   − −1 and n  (  ) ∼ N(0, Q  (  )) is the transition noise of the receiver state vector with covariance (  ) where  2   defines the continuous-time process noise intensity that has to be set based on the application with physical dimension [m 2 /s 3 ] and  2   the variance of the clock bias.This transition model is similar to the transition model presented in [29,30], with the advantage that the transition model is linear if the heading changes are known.
Since we incorporate only the heading changes Ψ(  ), we do not have any speed measurements and the speed has to be estimated implicitly.In order to adapt quickly to different walking speeds,    has to be large or many particles have to be used to cover all possible movements.A large value for    may cause backward movements of the transition model which results in estimation errors of Channel-SLAM.Another drawback of this transition model is that the estimated state information is completely lost when the user is standing; thus, the velocity components are zero.In order to overcome the mentioned problems, we develop a second transition model as described in the following.

Rician-Transition-Model.
Similarly to Section 3.2.1,we follow a two-dimensional positioning approach in the Cartesian coordinate system with the receiver state vector x  (  ) = [ ȓ  (  ), k  (  ),   (  )]  which consists of the - receiver positions ȓ  (  ) as defined in ( 18), the receiver velocity vector k  (  ), and the receiver's clock bias   (  ).The receiver velocity vector is where V  (  ) is the receiver speed and Ψ  (  ) the heading of the receiver; see Figure 4.The heading Ψ  (  ) describes the walking direction of the pedestrian with respect to the Cartesian coordinate system.Hence, we can define the transition model with where the velocity follows a Rician distribution with with scale parameter   (  ).For speeds close to zero, the Rician distribution approximates a Rayleigh distribution; thus, the speed is always positive.Hence, it has the advantage of preventing the filter from converging to negative velocities, which are highly unlikely regarding a normal pedestrian walking behavior (Please note that the developed transition model does not include standing or walking backwards phases.This could be additionally considered by extracting more information from the IMU measurements).This is important for our approach, since ambiguities of Channel-SLAM could otherwise cause a movement in the wrong direction.On the other hand, for higher speeds, the distribution becomes approximately Gaussian which reflects empirical data of pedestrian walking speeds [43].
Finally, the heading of the user is defined by where Ψ(  ) is the heading change from the IMU after calibration with the heading noise  Ψ (  ) using a von Mises distribution.For the transition prior of the clock bias we use similar to Section 3.2.1 a standard clock bias model   (  ) =   ( −1 ) +   ( −1 ), where   ( −1 ) defines the transition noise with the variance  2   [12,48] (Please note that if transmitter and receiver oscillators provide different frequencies, a clock drift parameter has to be considered additionally).In the following we refer to this transition model as Rician-Transition-Model.

Rao-Blackwellized Particle Filter. As introduced in [31],
Channel-SLAM is derived based on Rao-Blackwellization where the state space of x(  ) is partitioned into subspaces.
x (1)  VT,0 (t k ) Figure 5: The algorithm is based on a superordinate particle filter (superPF) and subordinate particle filters (subPFs).Each particle  = 1 ⋅ ⋅ ⋅   of the superPF consists of (  ) subPFs.Hence, we use particle filters (PFs) to estimate the subspaces representing the VTs inside a PF.The reason to use a PF instead of a low complexity EKF is the high nonlinearity of the measurements in (16).As shown in Figure 5, the algorithm is based on a superordinate particle filter (superPF) and subordinate particle filters (subPFs): Each particle  = 1 ⋅ ⋅ ⋅   of the superPF with the state vector x ()   (  ) = [r ()   (  )  , k ()  (  )  ,  ()  (  )]  consists of (  ) subPFs.Each subPF is represented by the particles x (,) VT, (  ) with  = 1, . . .,  ,, (  ), where  ,, (  ) stands for the number of particles in the th subPF with  = 0, . . ., (  ) − 1, estimating x () VT, (  ).Using subPFs for each VT allows using different numbers of particles in each subPF and, furthermore, allows a dynamic adjustment of the number of particles for each subPF introduced in Section 3.4.
According to [31], the marginalized posterior filtered density p(x  (  ) | z( 1: ), Ψ( 1: )) of the superPF can be approximated by importance samples (see [31,39]) as where  () (  ) defines the weight for the th particle at time instant   with and the weight Contrarily to [31], resampling is performed at each time instant to prevent degeneration; hence, ( 27) and ( 28) do not depend on the weights  () ( −1 ) and  (,)  ( −1 ).Additionally, the derivations in [31] consider a regularized PF [39] where x (,) VT, (  ) is drawn after resampling from the Gaussian-Kernel (⋅).The Gaussian-Kernel (⋅) improves the robustness of Channel-SLAM to cope with small model mismatches in the measurements.=1 are distributed on a grid inside a circle around r ()   ( 0 ) with radius d ( 0 ) + Δ  such that with spacing Δ  (the number of grid points  ,, (  ) can be estimated by Gauss's circle problem).The additional propagation length is VT, ( 0 ) − r ()  ( 0 )‖, where we inherently assume   ( 0 ) = 0 for the initialization.Hence, the total number of particles can be calculated as For each particle  of the superPF, the receiver state x ()  ( −1 ) is propagated according to the transition model described in Section 3.2 indicated by the function f(x ()  ( −1 ), Ψ(  ), n  (  )) using the heading changes Ψ(  ) (cf.Line (5) in Algorithm 1).Afterwards, Channel-SLAM determines whether the number of tracked MPCs has changed.In case that new MPCs have been detected, new subPFs are added and initialized using (29) (cf.Line (7)  if New paths detected then (7) Initialize new subPFs; (8) if Tracking of paths lost then (9) Delete corresponding subPFs; (10) for  = 0 : (  ) − 1 do (11) for  = 1 :  anymore, the corresponding subPFs are removed (cf.Line (9) in Algorithm 1).Equivalent to [31], neither KEST nor Channel-SLAM considers retracking of previous MPCs.Hence, if the tracking of an MPC has been lost and is regained again, the corresponding VT is initialized without any prior information.According to (14), the state x VT, (  ) is time-invariant; hence, each subPF assigns the states of the VTs with x (,) VT, ( −1 ).Thereafter, the weights for the subPFs and superPF are calculated using ( 28) and (27).
Afterwards, the subPFs and superPF are resampled.The basic idea of the resampling method is to eliminate particles with low weights and reproduce particles with high weights.Algorithm 2 shows a pseudocode of the resampling algorithm of Channel-SLAM which is based on the systematic resampling algorithm [49].Similarly to Algorithm 1, the Channel-SLAM resampling algorithm consists of a resampling algorithm for the superPF which includes resampling algorithms for the subPFs.First, the estimated sampled cumulative distribution function (CDF) of the superPF is constructed, presented by a vector c  with length   and element [c  ]  with  = 1, . . .,   .According to the estimated sampled CDF of the superPF, the subPFs are eliminated or resampled.The particles of the superPF with index  are assigned to the resampled particle index ; see Algorithm 2, Line (10), for the assignment of the receiver state.Afterwards, the (, )th subPF is resampled with  = 0, . . ., (  ) − 1 using a systematic resampling algorithm, where c  represents the estimated sampled CDF of the (, )th subPF.
As mentioned before, whenever a new MPC is tracked, many particles are initialized to cover all possible VT positions in each subPF.For example if the th MPC has a delay of d ( 0 ) = 30 m, each th subPF would use  ,, ( 0 ) = 2821 particles according to (29) for  = 1, . . .,   with spacing Δ  = 1 m.However during the receiver movement many particles of the subPFs are resampled at the same grid point because the states of the VTs x VT, (  ) are time-invariant.In order to adapt the number of particles, we limit the number of resampled particles per grid point to   ; see Algorithm 2, Line (20).Evaluations in Section 5 show that the reduction of the number of particles does not influence the positioning accuracy but, however, leads to a gain on computational performance.Afterwards, the states of the VTs x (,) VT, (  ) are drawn using a Gaussian-Kernel (cf.Line (26) in Algorithm 2).As stated in [31], the Gaussian-Kernel has a low bandwidth which does not influence the grid based reduction method.

Posterior Cramér-Rao Lower Bound
The PCRLB can be calculated by the inverse of the posterior information matrix J(  ) and provides a lower bound of the variance of a Bayesian estimator [50] with The inequality in (32) means that the difference M(  ) − J(  ) −1 is a positive semidefinite matrix.For the performance evaluation of a filter like Channel-SLAM with the system equations of ( 10) and ( 11), the posterior information matrix can be calculated recursively according to [51], with where x( −1 ) ln p (x (  ) | x ( −1 ) , Ψ (  ))] , where ∇  stands for the first-order partial derivatives with respect to  and Δ   stands for the second-order partial derivatives with Δ   ≜ ∇  ∇   .In order to calculate the PCRLB, we use the transition model of Section 3.2.1.In case of non-Gaussian noise and nonlinearity as in Section 3.2.2, the expectation estimator in (34) has to be approximated by Monte Carlo simulations.To calculate the PCRLB, we combine the time-invariant transition model for the VTs x VT (  ) as introduced in ( 14) and the transition model of the receiver (35), with where n  (  ) ∼ N(0, Q(  )) is the transition noise with covariance matrix where Q  (  ) is defined in (21).
Under the condition of a known control input Ψ(  ), see [51], and using the linear transition model of Section 3.2.1 and Gaussian distributed transition noise of (36), we obtain for (34) where the matrix F(  ) is the snapshot based Fisher information matrix.Substituting (37) into (33), we obtain using the matrix inversion lemma because of the singularity of Q(  ).
The snapshot based Fisher information matrix F(  ) in ( 38) can be obtained by with the covariance matrix R(  ), and with the propagation lengths   (  ) of ( 16) for the th MPC.

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

Conclusions
In this paper, we extended the work on multipath assisted positioning, called Channel-SLAM.The new positioning method takes advantage of the multipath components (MPCs) instead of mitigating them.In the proposed approach, multipath signals are treated as signals from virtual transmitters (VTs), where the locations of these VTs are unknown.To improve the accuracy of Channel-SLAM, an inertial measurement unit (IMU) is used to obtain heading information of the moving receiver.Furthermore, we investigate a novel particle filter (PF) implementation which adapts the number of particles during runtime.Measurement

3. 1 .
Position Estimation.Figure 2 presents the available sensors together with the corresponding measurements.As shown on the left, we measure the sampled received signal y(  ) as stated in (3) where we assume that the transmitter continuously emits known wideband signals.Based on y(  ), the multipath parameters amplitude   (  ) and delay   (  ) =   (  )/ for each MPC are estimated and tracked by KEST.The estimated propagation path lengths d (  ) = τ (  ) ⋅  of all (  ) MPCs of KEST are used as measurements z (  ) = [ d0 (  ) , . . ., d(  )−1 (  )]

Figure 2 :
Figure 2: System model consisting of a terrestrial receiver and an IMU.

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

Figure 4 :
Figure 4: Illustration of the prediction model for the pedestrian.

Figure 6 :Figure 7 :
Figure 6: Measurement scenario with a fixed transmitter and a moving receiver (pedestrian).The pedestrian moves on the blue track for 155 s, in total 111 m.The starting position and end position are indicated by the green and magenta circles.The metalized doors of the hangar and the chain-link fences act as reflecting surfaces for the transmitted wireless signal.Hence, we can calculate the corresponding VT positions by mirroring the physical transmitter position on the reflecting surfaces.

Figure 8 : 4 Figure 9 :
Figure 8: Recorded CIRs versus the pedestrian moving time in seconds.

Figure 11 :
Figure 11: Total number of particles   (  ) versus the pedestrian moving time in seconds for Dynamic-Channel-SLAM-1 and RBPF-Channel-SLAM-1.
, | 2 represents the square of the Frobenius norm of A.