Next Article in Journal
Mapping Arctic Lake Ice Backscatter Anomalies Using Sentinel-1 Time Series on Google Earth Engine
Previous Article in Journal
Determination of Venus’ Interior Structure with EnVision
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Robust Multipath-Assisted SLAM with Unknown Process Noise and Clutter Intensity

School of Electronics and Information, Northwestern Polytechnical University, Xi’an 710072, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(9), 1625; https://doi.org/10.3390/rs13091625
Submission received: 27 February 2021 / Revised: 15 April 2021 / Accepted: 19 April 2021 / Published: 21 April 2021

Abstract

:
In multipath-assisted simultaneous localization and mapping (SLAM), the geometric association of specular multipath components based on radio signals with environmental features is used to simultaneously localize user equipment and map the environment. We must contend with two notable model parameter uncertainties in multipath-assisted SLAM: process noise and clutter intensity. Knowledge of these two parameters is critically important to multipath-assisted SLAM, the uncertainty of which will seriously affect the SLAM accuracy. Conventional multipath-assisted SLAM algorithms generally regard these model parameters as fixed and known, which cannot meet the challenges presented in complicated environments. We address this challenge by improving the belief propagation (BP)-based SLAM algorithm and proposing a robust multipath-assisted SLAM algorithm that can accommodate model mismatch in process noise and clutter intensity. Specifically, we describe the evolution of the process noise variance and clutter intensity via Markov chain models and integrate them into the factor graph representing the Bayesian model of the multipath-assisted SLAM. Then, the BP message passing algorithm is leveraged to calculate the marginal posterior distributions of the user equipment, environmental features and unknown model parameters to achieve the goals of simultaneous localization and mapping, as well as adaptively learning the process noise variance and clutter intensity. Finally, the simulation results demonstrate that the proposed approach is robust against the uncertainty of the process noise and clutter intensity and shows excellent performances in challenging indoor environments.

Graphical Abstract

1. Introduction

Simultaneous localization and mapping (SLAM) [1,2] aims to estimate the time-varying position of the mobile agent when it is moving in an unknown environment and build a map of the surroundings from measurements provided by sensors. SLAM is significant in many fields, such as robotics [3], autonomous driving [4], and situation awareness [5]. In multipath-assisted methods [6,7,8,9,10], the relation between multipath components (MPCs) of radio signals and local geometry is used to turn the multipath propagation from an error source to an advantage that can improve the positioning accuracy. Recently, multipath-assisted SLAM, a combination of the two, has attracted considerable attention from scholars and has been applied in the field of indoor positioning [11,12,13,14,15,16].
Instead of mapping the environment in conventional SLAMs, the map in multipath-assisted SLAM is represented by an unknown number of features with unknown spatial positions. Features denote the physical transmitters (PTs) and virtual transmitters (VTs) in the environment. Radio signals are transmitted from base stations (PTs) to the user equipment (UE). Specular MPCs in the multipath radio channel are described by VTs, which are mirror images of the PTs. Due to the introduction of VTs, specular MPCs can be interpreted as line-of-sight (LOS) signals from the VTs, which can increase the number of transmitters involved in positioning and even realize localization in the case of insufficient PTs. Multipath-assisted SLAM detects the VTs associated with the PTs and estimates their positions as well as the trajectory of the UE. Multipath-assisted SLAM can cope with harsh multipath channel conditions and realize positioning and mapping without knowing the prior position of the PTs or without having a sufficient number of PTs, which greatly improves the accuracy and robustness of indoor positioning [17,18,19,20,21].
The data association (DA) between measurements and PTs or VTs is unknown, and measurements will be affected by noise, false alarms and missed detections, which makes the multipath-assisted SLAM complex. A multipath-assisted SLAM algorithm called Channel-SLAM was proposed in [22]. Channel-SLAM is based on a recursive Bayesian filter and can realize SLAM only with the prior position of the PT and the initial position and orientation of the UE. However, Channel-SLAM does not consider the DA uncertainty. Channel-SLAM is extended in [23] by adding a multiple hypothesis tracking (MHT)-based DA algorithm, which enables Channel-SLAM to associate the reappearing MPCs with corresponding VTs. Simulation results demonstrate that the improved Channel-SLAM algorithm has higher positioning reliability and lower computational complexity than the original method. A multipath-assisted SLAM algorithm called BP-SLAM based on belief propagation (BP) [24] was proposed in [25]. The approach jointly realizes the probabilistic DA and the estimation of UE positions and feature positions in the environment map. Simulation results demonstrate that BP-SLAM has lower complexity and higher scalability than other state-of-the-art SLAM algorithms and has higher positioning and mapping accuracy in the case of measurements containing false alarms and missed detections.
Furthermore, unknown and time-varying model parameters also complicate multipath-assisted SLAM. For example, the detection probability describes the possibility of PTs or VTs being detected by the UE. In practice, due to the influence of electromagnetic interference and surrounding environments, the detection probability is usually variable. However, the influence of the detection probability was not considered or simply considered to be fixed and known in [22,23,24,25]. Adaptive algorithms of the time-varying detection probability-based BP-SLAM algorithm were proposed in [26,27,28]. In [26], the detection probability was considered to be a discrete random variable with the value taking from a finite set. Then, the evolution of the detection probability was integrated into the multipath-assisted SLAM Bayesian model as an independent variable node. Thus, the algorithm can realize the online adaptive adjustment of the detection probability. In [27,28], the detection probability was expressed as a function of the MPC amplitude, and improved BP-SLAM algorithms were proposed, which leveraged MPC amplitude information to adaptively adjust the detection probability.
In multipath-assisted SLAM, we have to address the uncertainty of two important model parameters, process noise and clutter intensity, in addition to the detection probability. The change in process noise is mainly reflected in the change in UE dynamic models. Maneuvering or sudden acceleration and turning may cause process model uncertainties in UEs. Clutter represents spurious measurements that do not come from any transmitter (PTs or VTs). Clutter intensity is used to describe the amount of clutter in the environment. To the best of the authors’ knowledge, none of the multipath-assisted SLAM algorithms consider the uncertainty of process noise and clutter intensity at present. However, in the field of multitarget tracking (MTT) [29], which is closely related to feature-based SLAM, some algorithms are proposed with the condition that these two model parameters are unknown. In [30,31], clutter was regarded as measurements generated by incorrect targets. By distinguishing between correct and incorrect targets, the probability hypothesis density (PHD) filter, the cardinalized PHD (CPHD) filter, and the multi-Bernoulli filter were developed to independently estimate the clutter intensity in the environment. The conventional PHD filter and CPHD filter were extended in [32] to include estimation of the time-varying target birth intensity. In [33], an adaptive MTT algorithm based on the factor graph and the sum product algorithm [34] was proposed. The algorithm can be applied to complex scenes where the target dynamic model and detection probability change at the same time.
Most of the existing multipath-assisted SLAM algorithms regard the process noise and clutter intensity as known and fixed. However, in practical applications, the dynamic model of UEs cannot remain invariable all the time, and the process noise of UEs may change at any time. Clutter environments are often complicated due to the influence of electromagnetic interference, signal scattering and diffraction; hence, it is difficult to match the assumption of known process noise and clutter intensity with reality. Significant mismatches in the process and clutter model parameters seriously affect the SLAM accuracy, and even the estimation results diverge. Therefore, it is of critical importance to study the robust multipath-assisted SLAM algorithm when the process and clutter model parameters are unknown. Based on the framework of the BP-SLAM algorithm, in this paper, we propose a robust multipath-assisted SLAM algorithm, which can address the uncertainties of process and clutter model parameters in multipath-assisted SLAM. First, the UE state is extended to an augmented state including the process noise, and the evolution of the augmented state is described by new prediction and update rules compared to the conventional BP-SLAM algorithm. Second, we add a new variable to represent the clutter intensity in the Bayesian model of BP-SLAM and describe the evolution of the clutter intensity via a Markov chain model. Finally, we integrate the augmented state of the UE and the new variable node into a factor graph representing the Bayesian model of BP-SLAM and use the BP message passing algorithm to calculate marginal posterior distributions of the process noise, the clutter intensity, the UE and other features to realize the online adjustment of the process noise and the clutter intensity in multipath-assisted SLAM. The simulation results demonstrate that the proposed algorithm can adaptively learn the process noise and clutter intensity for SLAM, improving the robustness of the multipath-assisted SLAM system.
The remainder of this paper is organized as follows: Section 2 provides the environment map and the signal model. Section 3 describes the system model with unknown process noise and clutter intensity. Moreover, our robust multipath-assisted SLAM algorithm is proposed in Section 4. Then, the simulation results are reported in Section 5. Finally, Section 6 contains the conclusion of the paper.

2. Environment Map and Signal Model

Multipath-assisted SLAM combines the geometric parameters (delays, angles of arrival (AOAs) and angles of departure (AODs)) of the LOS signal in received radio signals with the geometric parameters of specular MPCs to localize the UE and build the environment map. These parameters are obtained according to the relationship between the positions of the UE and the PTs or VTs. The VT positions are mirror images of the PT positions with respect to the reflection surface. From the UE side, the reflected and virtual propagation paths originating from the VTs have the same AOA and delay; hence, the reflected path can be described as a LOS path between the VT and the UE. Regardless of how the UE moves, the VTs will remain static as long as the PTs are static. As an example, Figure 1 depicts an indoor environment map containing two PTs and some of the corresponding VTs. The red solid circle and the blue solid box denote PT1 and PT2 at a 1 ( 1 ) and a 1 ( 2 ) , respectively. The gray lines represent the doors, windows and walls in the room. They all have smooth surfaces that can reflect the radio signals from the PTs. The red hollow circles and blue hollow boxes represent the corresponding VTs, which are mirror images of PT1 and PT2. The green dot is the estimated position of the UE at u n , while the magenta line denotes the real trajectory of the UE.
We consider an environment map with one UE and J PTs. The UE position is unknown and time-varying, which is denoted by u n 2 , while the PT positions a 1 ( j ) 2 , j { 1 , , J } may be unknown, where the number of PT J is supposed to be known. Assume that there are L n ( j ) 1 VTs associated with the jth PT and that their positions are unknown and represented by a l ( j ) 2 , l { 2 , , L n ( j ) } . PTs and VTs are features of the environmental map, and their number is unknown and varies with the UE position. At each discrete time n, the radio signal s ( t ) is transmitted from the jth PT to the mobile UE, and the received signal of the UE can be expressed as [17]
s RX , n ( t ) = l = 1 L n ( j ) l , n ( j ) s ( t τ l , n ( j ) ) + g n ( j ) ( t ) + ω GWN ( t )
The first term on the right in Equation (1) describes L n ( j ) specular MPCs with complex amplitude l , n ( j ) and time delay τ l , n ( j ) , where l { 1 , , L n ( j ) } . These MPCs correspond to the features (PTs and VTs) of the environment. The delay τ l , n ( j ) corresponds to the distance between the UE and the jth PT ( l = 1 ) or the associated VTs ( l { 2 , , L n ( j ) } ); thus, we have τ l , n ( j ) = u n a l ( j ) / c , where c is the speed of light. Furthermore, g n ( j ) ( t ) in (1) denotes the diffuse MPCs, which is regarded as interference. Finally, ω GWN ( t ) represents additive white Gaussian noise.
At each time n, the radio channel estimator [35] processes the signal s RX , n ( t ) received by the UE to obtain M n ( j ) delay estimations of MPCs τ ^ m , n ( j ) , m n ( j ) = { 1 , , M n ( j ) } . The relation between the number of MPC delay estimations M n ( j ) and the number of specular MPCs L n ( j ) is as follows: on the one hand, some specular MPCs are unlikely to be detected by the radio channel estimator without generating delay estimations; on the other hand, some of the delay estimations cannot correspond to any specular MPCs. Therefore, there is no exact equality between M n ( j ) and L n ( j ) , and M n ( j ) mainly depends on the UE position u n and the environment. The estimated distances z m , n ( j ) , m { 1 , , M n ( j ) } of the mth MPC of PT j are obtained by multiplying the delay estimation τ ^ m , n ( j ) with the speed of light c, which are the measurements of the proposed algorithm used in this paper. For convenience, we define the stacked vectors of z m , n ( j ) , z n ( j ) [ z 1 , n ( j ) z M n ( j ) , n ( j ) ] T , z n [ z n ( 1 ) T z n ( J ) T ] T and z [ z 1 T z n T ] T and the stacked vectors of M n ( j ) , m n [ M n ( 1 ) M n ( J ) ] T and m = [ m 1 T m n T ] T .

3. System Model with Unknown Process Noise and Clutter Intensity

3.1. UE State and PF States

The UE state at time n is x ˜ n [ u n T v n T ] T and v n denotes the UE velocity vector. Potential features (PFs) are used to describe the features in the environment map (PTs and VTs) because the number of VTs is unknown and time-varying. PFs are indexed via (j, k), where j { 1 , , J } and k K n ( j ) { 1 , , K n ( j ) } (which implies that there are K n ( j ) PFs for each PT j). Note that the set of PFs includes the PT itself, and the number of PFs is unknown and varies with the UE position. We use the indicator variable e k , n ( j ) { 0 , 1 } to indicate whether PFs truly exist, where e k , n ( j ) = 1 means that PF (j, k) exists at time n, and vice versa. The indicator variable e k , n ( j ) is used to estimate the existence probability of PFs, which will be introduced in Section 4.1. The state of the PF (j, k) at time n is represented as an augmented state including the indicator variable e k , n ( j ) , i.e., y k , n ( j ) = [ a k , n ( j ) T e k , n ( j ) ] T , where a k , n ( j ) denotes the PF (j, k) position. We also define the stacked vectors of y k , n ( j ) , i.e., y n ( j ) = [ y 1 , n ( j ) T y K n ( j ) , n ( j ) T ] T , y n = [ y n ( 1 ) T y n ( J ) T ] T and y = [ y 1 T y n T ] T .
At any time n, two kinds of PFs exist: one is the legacy PF, which was established in the past, and the other is the new PF, which is established for the first time at time n. The augmented states of legacy PFs and new PFs for PT j are expressed as y _ k , n ( j ) [ a _ k , n ( j ) T e _ k , n ( j ) ] T , k K n - 1 ( j ) and y ¯ m , n ( j ) = [ a ¯ m , n ( j ) T e ¯ m , n ( j ) ] T , m n ( j ) , respectively. The number of new PFs is equal to the number of measurements M n ( j ) at time n; thus, the set and number of legacy PFs are updated as
K n ( j ) = K n - 1 ( j ) n ( j ) , K n ( j ) = K n - 1 ( j ) + M n ( j )
From Equation (2), we see that the PF set updates as a dynamic increasing process, i.e., legacy PFs at time n − 1 and new PFs at time n constitute the legacy PFs at time n + 1. (Note that the number of PFs K n ( j ) at time n does not actually increase M n ( j ) compared with that at time n − 1, since some PFs with low existence probability are pruned, which will be explained in Section 4.1.) We also define the stacked vectors related to PF states as follows: for legacy PFs for PT j, a _ n ( j ) [ a _ 1 , n ( j ) T a _ K n 1 ( j ) , n ( j ) T ] T , e _ n ( j ) [ e _ 1 , n ( j ) e _ K n 1 ( j ) , n ( j ) ] T and y _ n ( j ) [ y _ 1 , n ( j ) T y _ K n 1 ( j ) , n ( j ) T ] T ; for new PFs for PT j, a ¯ n ( j ) [ a ¯ 1 , n ( j ) T a ¯ M n ( j ) , n ( j ) T ] T , e ¯ n ( j ) [ e ¯ 1 , n ( j ) e ¯ M n ( j ) , n ( j ) ] T and y ¯ n ( j ) [ y ¯ 1 , n ( j ) T y ¯ M n ( j ) , n ( j ) T ] T ; and for the combination of legacy PFs and new PFs for PT j, y n ( j ) [ y _ n ( j ) T y ¯ n ( j ) T ] T . Note that the entries of y n ( j ) are y k , n ( j ) and k K n ( j ) ; for all legacy PFs and new PFs, y _ n [ y _ n ( 1 ) T y _ n ( J ) T ] T and y ¯ n [ y ¯ n ( 1 ) T y ¯ n ( J ) T ] T .
For each PT, the measurements z m , n ( j ) are affected by the DA (measurement source) uncertainty. In other words, it is hard to know which measurement z m , n ( j ) corresponds to which PF k K n ( j ) . The measurement z m , n ( j ) may not originate from any PF (false alarms or clutter), or it may not be detected by the channel estimator (missed detections). We assume that each PF generates at most one measurement at any time n and that each measurement can be generated by at most one PF [36]. At time n, the association between the measurements z m , n ( j ) and PFs can be expressed as a feature-oriented K n - 1 ( j ) -dimensional DA vector h n ( j ) = [ h 1 , n ( j ) h K n - 1 ( j ) , n ( j ) ] T . If the legacy PF (j, k) generates a measurement z m , n ( j ) , the kth entry of h n ( j ) is h k , n ( j ) = m n ( j ) ; if it does not generate any measurement, then we have h k , n ( j ) = 0 . Furthermore, we consider a measurement-oriented M n ( j ) -dimensional DA vector q n ( j ) = [ q 1 , n ( j ) q M n ( j ) , n ( j ) ] T . If the measurement z m , n ( j ) originates from the legacy PF (j, k), the mth entry of q n ( j ) is q m , n ( j ) = k K n ( j ) ; if the measurement z m , n ( j ) does not originate from any legacy PF, then q k , n ( j ) = 0 . Additionally, we define the stacked vectors of DA vectors h n [ h n ( 1 ) T h n ( J ) T ] T , h [ h 1 T h n T ] T , q n [ q n ( 1 ) T q n ( J ) T ] T and q [ q 1 T q n T ] T . The DA vectors h and q are random and completely equivalent since one of them can be determined by the other.

3.2. Markov Chain Modeling of the Process Noise and the Clutter Intensity

In addition to the essential information to be tracked, i.e., the UE state, PF states and DA vectors, model parameters such as the process noise and the clutter intensity are also unknown and time-varying. The process noise determines the difference of the UE dynamic models (assuming that the UE is in a nearly constant velocity model), and the clutter intensity describes the amount of clutter (false alarms) in the environment. Significant mismatches in the process and clutter model will seriously affect the SLAM accuracy. To obtain a more robust multipath-assisted SLAM algorithm, we propose tracking unknown model parameters, process noise and clutter intensity in the framework of the BP-SLAM algorithm. The process noise and the clutter intensity can be modeled in the same way: ρ n ( s ) , s = 1 , , S is assumed to be an unknown model parameter. We model ρ n ( s ) as a time-varying discrete random variable taking values from a finite set A = { ω 1 ( s ) , , ω N s ( s ) } and assume that the distribution of the initial parameter ρ 1 ( s ) corresponds to a probability mass function (PMF) p ( ρ 1 ( s ) ) , ρ 1 ( s ) A . Additionally, we assume that the evolution ρ n ( s ) is independent, which corresponds to the first-order Markov chain model with transition matrices P s [ 0 , 1 ] N s × N s , where [ 0 , 1 ] N s × N s refers to the elements in the matrix with size N s × N s being [ 0 , 1 ] . Thus, the PMF of the state transition of ρ n ( s ) is
p ( ρ n ( s ) = ω i ( s ) | ρ n 1 ( s ) = ω i ( s ) ) = [ P s ] i , i
Note that i = 1 N s [ P s ] i , i = 1 for all i = 1 , , N s . We define the stacked vector of ρ n ( s ) ρ [ ρ 1 T ρ n T ] T , where ρ n [ ρ n ( 1 ) ρ n ( S ) ] T . Based on the assumptions above, the prior PMF of ρ factorizes as
p ( ρ ) = s = 1 S p ( ρ 1 ( s ) ) n = 2 n p ( ρ n ( s ) | ρ n 1 ( s ) )

3.3. State Evolution with Unknown Process Noise

The probability density function (PDF) of the UE state transition is determined by the process noise variance if the UE is in a nearly constant velocity model. Using the interacting multiple model (IMM) algorithm [37], the UE can switch among nearly constant velocity models with different process noise variances at any time. To track the change of the UE dynamic model, the original state of the UE x ˜ n is augmented by the UE process noise variance as x n [ x ˜ n T σ n 2 ] T , where σ n 2 denotes the process noise variance of the UE at time n. We model the process noise variance σ n 2 as a discrete random variable taking values from a finite set = { θ 1 , , θ N b } , the evolution of which is based on the Markov chain model of Section 3.2 with the transition matrix Q b [ 0 , 1 ] N b × N b . Thus, the PMF of the state transition of σ n 2 is
p ( σ n 2 = θ i | σ n 1 2 = θ i ) = [ Q b ] i , i , θ i , θ i
Additionally, we define the stacked vectors σ [ σ 1 2 σ n 2 ] T and x [ x 1 T x n T ] . The evolutions of the UE augmented state x n and the legacy PF augmented state y _ k , n ( j ) are assumed to be independent and based on the Markov state model. Their evolutions can be jointly expressed as
f ( x n , y _ n | x n 1 , y n 1 ) = f ( x ˜ n , σ n 2 | x ˜ n - 1 , σ n 1 2 ) f ( y _ n | y n 1 ) = f ( x ˜ n | σ n 2 , x ˜ n - 1 ) p ( σ n 2 | σ n 1 2 ) × j = 1 J k = 1 K n 1 ( j ) f ( y _ k , n ( j ) | y k , n 1 ( j ) )
where f ( x ˜ n | σ n 2 , x ˜ n 1 ) denotes the state transition PDF of the UE with process noise variance σ n 2 , and we assume that the UE original states x ˜ n and σ n 2 are statistically independent across time n here; f ( y _ k , n ( j ) | y k , n 1 ( j ) ) = f ( a _ k , n ( j ) , e _ k , n ( j ) | a k , n 1 ( j ) , e k , n 1 ( j ) ) denotes the state transition PDF of the legacy PF (j, k). If PF (j, k) exists at time n−1 ( e k , n 1 ( j ) = 1 ), it either disappears ( e _ k , n ( j ) = 0 ) or continues to exist ( e _ k , n ( j ) = 1 ) at time n. In the latter case, PF (j, k) becomes a legacy PF at time n. If PF (j, k) survives at time n, its update can be represented by the state transition PDF f ( a _ k , n ( j ) | a k , n 1 ( j ) ) . Thus, when e k , n 1 ( j ) = 1 , f ( y _ k , n ( j ) | y k , n 1 ( j ) ) in Equation (6) is expressed as
f ( a _ k , n ( j ) , e _ k , n ( j ) | a k , n 1 ( j ) , e k , n 1 ( j ) = 1 ) = { ( 1 P s ) f D ( a _ k , n ( j ) ) , e _ k , n ( j ) = 0 P s f ( a _ k , n ( j ) | a k , n 1 ( j ) ) , e _ k , n ( j ) = 1
where P s denotes the survival probability of the PF and f D ( a _ k , n ( j ) ) denotes an arbitrary “dummy PDF” [38]. If PF (j, k) does not exist at time n − 1, it cannot exist at time n. Thus, we have
f ( a _ k , n ( j ) , e _ k , n ( j ) | a k , n 1 ( j ) , e k , n 1 ( j ) = 0 ) = { f D ( a _ k , n ( j ) ) , e _ k , n ( j ) = 0 0 , e _ k , n ( j ) = 1

3.4. Prior Distributions with Unknown Clutter Intensity

Clutter are false measurements in the environment, which are usually modeled as Poisson point processes; hence, all the characteristics of the clutter can be completely described by the intensity information. The clutter intensity is defined as the product of the average number of clutter at each time λ CL , n ( j ) and the spatial distribution PDF of clutter f CL ( z m , n ( j ) ) . f CL ( z m , n ( j ) ) is generally considered to be a priori; therefore, the clutter intensity can only be determined by the average number of clutter λ CL , n ( j ) . λ CL , n ( j ) is assumed to be statistically independent across j and modeled as a discrete random variable with values taken from a finite set G = { r 1 , , r N g } . The evolution of λ CL , n ( j ) is based on the Markov chain model of Section 3.2 with the transition matrix G ( j ) [ 0 , 1 ] N g × N g , where G ( j ) is equal at all times and only related to PT. The initial distribution of λ CL , n ( j ) is given by the PMF p ( λ CL , 1 ( j ) ) . We also define the stacked vector λ CL [ λ CL , 1 T λ CL , n T ] T , where λ CL , n [ λ CL , 1 ( j ) λ CL , n ( j ) ] T . According to (2), the prior PMF of λ CL factorizes as
p ( λ CL ) = j = 1 J p ( λ CL , 1 ( j ) ) n = 2 n p ( λ CL , n ( j ) | λ CL , n 1 ( j ) )
where p ( λ CL , n ( j ) = r i | λ CL , n 1 ( j ) = r i ) = [ G ( j ) ] i , i . Similar to clutter, we model the number of newly detected features as a Poisson point process with mean value λ new , n ( j ) . We also define the detection probability P d ( j ) ( x n , a k , n ( j ) ) , which refers to the probability that the PF can be detected or the probability that the PF can generate the measurement z m , n ( j ) in the channel estimation process [39]. The detection probability P d ( j ) ( x n , a k , n ( j ) ) ( 0 , 1 ] is supposed known in this paper.
Then, given the UE augmented state x n , the legacy PF augmented state y _ n ( j ) and the average number of clutter λ CL , n ( j ) , the joint prior conditional PMF of the association vectors h n ( j ) and q n ( j ) , the existence indicator variable of new PFs e ¯ n ( j ) and the number of measurements M n ( j ) is expressed as [38]
p ( h n ( j ) , q n ( j ) , e ¯ n ( j ) , M n ( j ) | x n , y _ n ( j ) , λ CL , n ( j ) ) ψ ( h n ( j ) , q n ( j ) ) ( m = 1 M n ( j ) c 1 ( e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; K n ( j ) ) ) × ( k = 1 K n 1 ( j ) b 1 ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) ; M n ( j ) ) )
The exclude indicator function in (10) is
ψ ( h n ( j ) , q n ( j ) ) k = 1 K n ( j ) m = 1 M n ( j ) φ ( h k , n ( j ) , q m , n ( j ) )
where when h k , n ( j ) = m , q m , n ( j ) k or q m , n ( j ) = k , h k , n ( j ) m , φ ( h k , n ( j ) , q m , n ( j ) ) is 0; otherwise, φ ( h k , n ( j ) , q m , n ( j ) ) is 1. The indicator function enforces the assumption in Section 3.1 that each PF generates at most one measurement at any time n and that each measurement can be generated by at most one PF. Note that this DA constraint function is redundant since it considers both the feature-oriented and measurement-oriented association vectors h n ( j ) and q n ( j ) . However, such a redundant representation is the key to obtaining an algorithm that has low computational complexity and excellent scalability [38,40]. Furthermore, c 1 ( e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; K n ( j ) ) and b 1 ( x n , a _ k , n ( j ) , e _ k , n ( j ) , p k , n ( j ) ; M n ( j ) ) in (10) are defined as
c 1 ( 1 , q m , n ( j ) , λ CL , n ( j ) ; K n ( j ) ) { 0 , q m , n ( j ) K n ( j ) λ new , n ( j ) λ CL , n ( j ) , q m , n ( j ) = 0
c 1 ( 0 , q m , n ( j ) , λ CL , n ( j ) ; K n ( j ) ) f D ( a ¯ m , n ( j ) )
b 1 ( x n , a _ k , n ( j ) , 1 , h k , n ( j ) , λ CL , n ( j ) ; M n ( j ) ) { P d ( j ) ( x n , a _ k , n ( j ) ) λ CL , n ( j ) , h k , n ( j ) n ( j ) 1 P d ( j ) ( x n , a _ k , n ( j ) ) , h k , n ( j ) = 0
b 1 ( x n , a _ k , n ( j ) , 0 , h k , n ( j ) , λ CL , n ( j ) ; M n ( j ) ) 1 ( h k , n ( j ) )
where the indicator function 1 ( h k , n ( j ) ) is defined as 1 ( h k , n ( j ) ) = 1 when h k , n ( j ) = 0 ; otherwise, 1 ( h k , n ( j ) ) = 0 . In particular, when n = 1, since the legacy PF y _ 1 ( j ) is an empty vector, we have p ( h 1 ( j ) , q 1 ( j ) , e ¯ 1 ( j ) , M 1 ( j ) | x 1 , y _ 1 ( j ) , λ CL , 1 ( j ) ) = p ( h 1 ( j ) , q 1 ( j ) , e ¯ 1 ( j ) , M 1 ( j ) | x 1 , λ CL , 1 ( j ) ) .
Assume that the spatial distribution PDF of newly detected features is f new , n ( a ¯ m , n ( j ) ) (the calculation of which will be discussed in Section 5.1); then, for PT j, given e ¯ n ( j ) and M n ( j ) , the prior PDF of the new PF state a ¯ n ( j ) is
f ( a ¯ n ( j ) | e ¯ n ( j ) , M n ( j ) ) = ( m N e ¯ n ( j ) f new , n ( a ¯ m , n ( j ) ) ) m N ¯ e ¯ n ( j ) f D ( a ¯ m , n ( j ) )
where N ¯ e ¯ n ( j ) n ( j ) \ N e ¯ n ( j ) and N e ¯ n ( j ) { m n ( j ) : e ¯ m , n ( j ) = 1 } .
Suppose that the initial position of the UE is a priori, i.e., f ( x 1 ) is known. Then, the joint prior PDF of the UE augmented state x , the PF augmented state y , DA vectors h and q , the number of measurements m and the average number of clutters λ CL is expressed as
f ( x , y , h , q , m , λ CL ) = f ( x 1 ) ( j = 1 J p ( λ CL , 1 ( j ) ) f ( y ¯ 1 ( j ) , h 1 ( j ) , q 1 ( j ) , M 1 ( j ) | x 1 , λ CL , 1 ( j ) ) ) × n = 2 n f ( x ˜ n | σ n 2 , x ˜ n - 1 ) p ( σ n 2 | σ n 1 2 ) × j = 1 J ( k = 1 K n 1 ( j ) f ( y _ k , n ( j ) | y k , n 1 ( j ) ) ) p ( λ CL , n ( j ) | λ CL , n 1 ( j ) ) × f ( a ¯ n ( j ) | e ¯ n ( j ) , M n ( j ) ) p ( h n ( j ) , q n ( j ) , e ¯ n ( j ) , M n ( j ) | x n , y _ n ( j ) , λ CL , n ( j ) )
where f ( y ¯ 1 ( j ) , h 1 ( j ) , q 1 ( j ) , M 1 ( j ) | x 1 , λ CL , 1 ( j ) ) = p ( h 1 ( j ) , q 1 ( j ) , e ¯ 1 ( j ) , M 1 ( j ) | x 1 , λ CL , 1 ( j ) ) f ( a ¯ 1 ( j ) | e ¯ 1 ( j ) , M 1 ( j ) ) (see Equation (10) and Equation (16)), and the last four factors are obtained from Equation (6), Equation (9), Equation (10) and Equation (16), respectively.

3.5. Measurement Model and Likelihood Function

The relation between the measurement z m , n ( j ) and the position-dependent states x ˜ n and a k , n ( j ) is described by the conditional PDF f ( z m , n ( j ) | x ˜ n , a k , n ( j ) ) , an example of which will be given in Section 5.1. This PDF is the central element of the likelihood function f ( z n ( j ) | x n , y _ n ( j ) , y ¯ n ( j ) , h n ( j ) , M n ( j ) ) , which is a function of x n , y _ n ( j ) , y ¯ n ( j ) , h n ( j ) and M n ( j ) when the measurement z n ( j ) is observed. Assuming that z n ( j ) is observed, i.e., that z n ( j ) is fixed, then, up to a constant factor, the likelihood function can be expressed as [25]
f ( z n ( j ) | x n , y _ n ( j ) , y ¯ n ( j ) , h n ( j ) , M n ( j ) ) m N e ¯ n ( j ) f ( z m , n ( j ) | x ˜ n , a ¯ m , n ( j ) ) f CL ( z m , n ( j ) ) × ( k = 1 K n ( j ) b 2 ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) ; z n ( j ) ) )
where b 2 ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) ; z n ( j ) ) is defined as
b 2 ( x n , a _ k , n ( j ) , 1 , h k , n ( j ) ; z n ( j ) ) { f ( z m , n ( j ) | x ˜ n , a _ k , n ( j ) ) f CL ( z m , n ( j ) ) , h k , n ( j ) = m 1 , h k , n ( j ) = 0
and b 2 ( x n , a _ k , n ( j ) , 0 , p k , n ( j ) ; z n ( j ) ) 1 . In particular, at the initial time when n = 1, we have f ( z 1 ( j ) | x 1 , y _ 1 ( j ) , y ¯ 1 ( j ) , h 1 ( j ) , M 1 ( j ) ) = f ( z 1 ( j ) | x 1 , y ¯ 1 ( j ) , h 1 ( j ) , M 1 ( j ) ) . Based on Equation (18), the likelihood function of the stacked vector z of the measurement z m , n ( j ) involving all PTs j = 1 , , J and all times n = 1 , , n is
f ( z | x , y , h , m ) j = 1 J ( m N e ¯ 1 ( j ) f ( z m , 1 ( j ) | x ˜ 1 , a ¯ m , 1 ( j ) ) f CL ( z m , 1 ( j ) ) ) × n = 2 n ( m N e ¯ n ( j ) f ( z m , n ( j ) | x ˜ n , a ¯ m , n ( j ) ) f CL ( z m , n ( j ) ) ) × k = 1 K n ( j ) b 2 ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) ; z n ( j ) )
Now, we consider the product of the likelihood function Equation (18), the prior PDF of the association vectors Equation (10) and the prior PDF of new PFs Equation (16), and we have
f ( z n ( j ) | x n , y _ n ( j ) , y ¯ n ( j ) , h n ( j ) , M n ( j ) ) f ( a ¯ n ( j ) | e ¯ n ( j ) , M n ( j ) ) × p ( h n ( j ) , q n ( j ) , e ¯ n ( j ) , M n ( j ) | x n , y _ n ( j ) , λ CL , n ( j ) ) ψ ( h n ( j ) , q n ( j ) ) ( k = 1 K n 1 ( j ) b ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) ) × m = 1 M n ( j ) c ( x n , a ¯ m , n ( j ) , e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) )
where b ( , , , , ; z n ( j ) ) b 1 ( , , , , ; M n ( j ) ) b 2 ( , , , ; z n ( j ) ) and c ( , , , , ; z n ( j ) ) c 1 ( , , ; K n ( j ) ) c 2 ( , , , ; z n ( j ) ) . b 1 ( , , , , ; M n ( j ) ) , b 2 ( , , , ; z n ( j ) ) and c 1 ( , , ; K n ( j ) ) are given by Equation (14), Equation (12) and Equation (19), respectively, and c 2 ( x n , a ¯ m , n ( j ) , e _ m , n ( j ) , q m , n ( j ) ; z n ( j ) ) is defined as
c 2 ( x n , a ¯ m , n ( j ) , 1 , q m , n ( j ) ; z n ( j ) ) { 0 , q m , n ( j ) = k f new , n ( a ¯ m , n ( j ) ) f ( z m , n ( j ) | x ˜ n , a ¯ m , n ( j ) ) f CL ( z m , n ( j ) ) , q m , n ( j ) = 0
and c 2 ( x n , a ¯ m , n ( j ) , 0 , q m , n ( j ) ; z n ( j ) ) f D ( a ¯ m , n ( j ) ) . Then, b ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) and c ( x n , a ¯ m , n ( j ) , e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) can be expressed as
b ( x n , a _ k , n ( j ) , 1 , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) = { P d ( j ) ( x n , a _ k , n ( j ) ) f ( z m , n ( j ) | x ˜ n , a _ k , n ( j ) ) λ CL , n ( j ) f CL ( z m , n ( j ) ) , h k , n ( j ) = m 1 P d ( j ) ( x n , a _ k , n ( j ) ) , h k , n ( j ) = 0
and b ( x n , a _ k , n ( j ) , 0 , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) = 1 ( h k , n ( j ) ) ;
c ( x n , a ¯ m , n ( j ) , 1 , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) = { 0 , q m , n ( j ) = k λ new , n ( j ) f new , n ( a ¯ m , n ( j ) ) f ( z m , n ( j ) | x ˜ n , a ¯ m , n ( j ) ) λ CL , n ( j ) f CL ( z m , n ( j ) ) , q m , n ( j ) = 0
and c ( x n , a ¯ m , n ( j ) , 0 , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) = f D ( a ¯ m , n ( j ) ) .

4. The Proposed Algorithm

In this section, the proposed BP-based robust multipath-assisted SLAM algorithm is described in detail.

4.1. Detection and Estimation

The ultimate goal of our algorithm is to estimate the UE original state x ˜ n and PF states a k , n ( j ) from the measurements z involving all PTs and all times with unknown process noise and clutter intensity. In the framework of the Bayesian theorem, estimations of the UE and PF states are based on their respective posterior PDFs. First, to estimate x ˜ n , we use the minimum mean square error (MMSE) estimator, i.e.,
x ˜ ^ n MMSE = x ˜ n f ( x ˜ n | z ) d x ˜ n
To obtain f ( x ˜ n | z ) , we have to marginalize the posterior PDF of the UE augmented state f ( x n | z ) = f ( x ˜ n , σ n 2 | z ) , i.e.,
f ( x ˜ n | z ) = σ n 2 f ( x ˜ n , σ n 2 | z )
Furthermore, detecting whether PFs exist at time n depends on the posterior existence probability p ( e k , n ( j ) = 1 | z ) , which is obtained by marginalizing the posterior PDF of PF augmented states
p ( e k , n ( j ) = 1 | z ) = f ( a k , n ( j ) , e k , n ( j ) = 1 | z ) d a k , n ( j )
Then, if p ( e k , n ( j ) = 1 | z ) > P th , PF k K n ( j ) is detected at time n, where P th is the detection threshold. Thus, the estimation of the state of the detected PF a k , n ( j ) is
a ^ k , n ( j ) MMSE = a k , n ( j ) f ( a k , n ( j ) | e k , n ( j ) = 1 , z ) d a k , n ( j )
where
f ( a k , n ( j ) | e k , n ( j ) = 1 , z ) = f ( a k , n ( j ) , e k , n ( j ) = 1 | z ) p ( e k , n ( j ) = 1 | z )
Therefore, to estimate x ˜ n and a k , n ( j ) , we must calculate the marginal posterior PDFs f ( x n | z ) and f ( y k , n ( j ) | z ) of the augmented states of the UE and PFs. To prevent the unlimited increase in the number of PFs caused by Equation (2) in Section 3.1, namely, K n ( j ) = K n - 1 ( j ) + M n ( j ) , the posterior existence probability p ( e k , n ( j ) = 1 | z ) is also used to prune some PFs. Specifically, PFs are preserved only when p ( e k , n ( j ) = 1 | z ) exceeds the predefined prune threshold P pr .

4.2. Joint Posterior Distribution and Factor Graph

The posterior PDFs f ( x n | z ) and f ( y k , n ( j ) | z ) related in Equations (25)–(29) are marginal distributions of the joint posterior PDF f ( x , y , h , q , m , λ CL | z ) , which involves the UE augmented state x , PFs augmented states y , DA vectors h and q , the number of measurements m and the average number of clutters λ CL at all times. By running the BP message passing algorithm on the factor graph representing the factorization of the joint posterior PDF f ( x , y , h , q , m , λ CL | z ) , effective estimations of the corresponding marginal posterior PDFs can be obtained. Based on the Bayesian theorem and related independence hypotheses, using Equation (6), Equation (17), Equation (20) and Equation (21) and some simple derivations, the joint posterior PDF factorizes as
f ( x , y , h , q , m , λ CL | z ) f ( z | x , y , h , m ) f ( x , y , h , q , m , λ CL ) = f ( x 1 ) ( j = 1 J p ( λ CL , 1 ( j ) ) m = 1 M 1 ( j ) c ( x 1 , a ¯ m , 1 ( j ) , e ¯ m , 1 ( j ) , q m , 1 ( j ) , λ CL , 1 ( j ) ; z 1 ( j ) ) ) × n = 2 n f ( x n | x n 1 ) j = 1 J p ( λ CL , n ( j ) | λ CL , n 1 ( j ) ) ψ ( h n ( j ) , q n ( j ) ) × ( k = 1 K n 1 ( j ) f ( y _ k , n ( j ) | y k , n 1 ( j ) ) b ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) ) × m = 1 M n ( j ) c ( x n , a ¯ m , n ( j ) , e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) )
This factorization is represented by the factor graph [34] in Figure 2, in which the variable node labeled x denotes the UE augmented state x n [ x ˜ n T σ n 2 ] T , and the average number of clutter λ CL , n ( j ) is denoted by the independent variable node labeled λ CL .

4.3. BP Message Passing Algorithm

Now, we calculate the marginal posterior PDFs f ( x n | z ) and f ( y k , n ( j ) | z ) from Section 4.1 by running the BP message passing algorithm in the factor graph in Figure 2. Since the factor graph has loops, the computation sequence of each message must be specified. The order of the message calculation obeys the following rules: (1) the message in the current time is passed only to the next time; (2) the iterative message passing is conducted only in the DA section and executed once for each PT at all times; and (3) if an edge connects the variable nodes of the UE augmented state and new PF augmented states, the message can only be passed from the former to the latter. For the calculation equations of the message passing algorithm, we need to consider the BP message when PFs do not exist ο ( y k , n ( j ) ) = ο ( a k , n ( j ) , e k , n ( j ) ) , e k , n ( j ) = 0 . Similar to the “dummy PDF” introduced in Section 3.3, we define ο ( a k , n ( j ) , 0 ) = ο k , n ( j ) (note that this message is a single number, not a PDF). Combining the calculation order of messages and the computation method of BP messages, the steps of the BP message passing algorithm are described as follows:
  • Prediction.
First, the prediction message of the UE augmented state α ( x n ) = α ( x ˜ n , σ n 2 ) is expressed as
α ( x ˜ n , σ n 2 ) = σ n 1 2 f ( x ˜ n , σ n 2 | x ˜ n 1 , σ n 1 2 ) f ( x ˜ n 1 , σ n 1 2 ) d x ˜ n 1
Substituting (6) into (31), we have
α ( x ˜ n , σ n 2 ) = σ n 1 2 f ( x ˜ n | σ n 2 , x ˜ n 1 ) p ( σ n 2 | σ n 1 2 ) f ( x ˜ n 1 , σ n 1 2 ) d x ˜ n 1
where f ( x ˜ n 1 , σ n 1 2 ) is the calculated belief of the UE augmented state at time n − 1. Then, the prediction message of legacy PF augmented states y _ k , n ( j ) [ a _ k , n ( j ) T e _ k , n ( j ) ] T is
α k ( a _ k , n ( j ) , e _ k , n ( j ) ) = e k , n - 1 ( j ) { 0 , 1 } f ( a _ k , n ( j ) , e _ k , n ( j ) | a k , n 1 ( j ) , e k , n 1 ( j ) ) f _ ( a k , n 1 ( j ) , e k , n 1 ( j ) ) d a k , n 1 ( j )
where the belief of the legacy PF augmented state f _ ( a k , n 1 ( j ) , e k , n 1 ( j ) ) is also calculated at time n − 1. Substituting (7) and (8) into (33), when e _ k , n ( j ) = 1 , we obtain
α k ( a _ k , n ( j ) , 1 ) = P s f ( a _ k , n ( j ) | a k , n 1 ( j ) ) f _ ( a k , n 1 ( j ) , 1 ) d a k , n 1 ( j )
while when e _ k , n ( j ) = 0 ,
α k , n ( j ) = ( 1 - P s ) f _ ( a k , n 1 ( j ) , 1 ) d a k , n 1 ( j ) + f _ k , n 1 ( j )
where α k , n ( j ) α k ( a _ k , n ( j ) , 0 ) d a _ k , n ( j ) and f _ k , n 1 ( j ) f _ ( a k , n 1 ( j ) , 0 ) d a k , n 1 ( j ) . Finally, the prediction also includes the prediction of the average number of clutter, which is calculated as
χ ( λ CL , n ( j ) ) = λ CL , n 1 ( j ) G p ( λ CL , n ( j ) | λ CL , n 1 ( j ) ) p ( λ CL , n 1 ( j ) )
Similarly, p ( λ CL , n 1 ( j ) ) denotes the belief of the average number of clutter at time n − 1. After the prediction step, the computation of other messages will be performed with all PTs j { 1 , , J } , all legacy PFs k K n 1 ( j ) and all new PFs m M n ( j ) in parallel.
2.
Measurement evaluation of legacy PFs.
In Figure 2, the message ε ( h k , n ( j ) ) passed to the variable node representing the feature-oriented DA vector h k , n ( j ) is calculated as
ε ( h k , n ( j ) ) = λ CL , n ( j ) G σ n 2 χ ( λ CL , n ( j ) ) α k ( a _ k , n ( j ) , 1 ) α ( x ˜ n , σ n 2 ) × b ( x n , a _ k , n ( j ) , 1 , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) d x ˜ n d a _ k , n ( j ) + 1 ( h k , n ( j ) ) α k , n ( j )
3.
Measurement evaluation of new PFs.
The calculation of the message ξ ( q m , n ( j ) ) passed to the variable node representing the measurement-oriented DA vector q m , n ( j ) is
ξ ( q m , n ( j ) ) = λ CL , n ( j ) G e ¯ m , n ( j ) { 0 , 1 } σ n 2 χ ( λ CL , n ( j ) ) α ( x ˜ n , σ n 2 ) × c ( x n , a ¯ m , n ( j ) , e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) d x ˜ n d a ¯ m , n ( j )
Substituting (24) into (38), ξ ( q m , n ( j ) ) can be expressed as
ξ ( q m , n ( j ) ) = { 1 , q m , n ( j ) K n 1 ( j ) 1 + λ CL , n ( j ) G σ n 2 λ new , n ( j ) λ CL , n ( j ) f CL ( z m , n ( j ) ) χ ( λ CL , n ( j ) ) × α ( x ˜ n , σ n 2 ) f new , n ( a ¯ m , n ( j ) ) f ( z m , n ( j ) | x ˜ n , a ¯ m , n ( j ) ) d x ˜ n d a ¯ m , n ( j ) , q m , n ( j ) = 0
4.
Iterative DA
Next, using ε ( h k , n ( j ) ) and ξ ( q m , n ( j ) ) , the messages γ ( h k , n ( j ) ) and ν ( q m , n ( j ) ) can be calculated by the iterative BP method. First, for each measurement m n ( j ) , the iterative calculations of the messages δ m k ( t ) ( h k , n ( j ) ) and ς k m ( t ) ( q m , n ( j ) ) are [41]
δ m k ( t ) ( h k , n ( j ) ) = q m , n ( j ) = 0 K n 1 ( j ) ξ ( q m , n ( j ) ) φ ( h k , n ( j ) , q m , n ( j ) ) k K n 1 ( j ) \ { k } ς k m ( t 1 ) ( q m , n ( j ) )
ς k m ( t ) ( q m , n ( j ) ) = h k , n ( j ) = 0 M n ( j ) ε ( h k , n ( j ) ) φ ( h k , n ( j ) , q m , n ( j ) ) m n ( j ) \ { m } δ m k ( t ) ( h k , n ( j ) )
where the number of iterations is t = 1 , , T . When i = 0 , the initial value of the iterations in (40) and (41) is
ς k m ( 0 ) ( q m , n ( j ) ) = h k , n ( j ) = 0 M n ( j ) ε ( h k , n ( j ) ) φ ( h k , n ( j ) , q m , n ( j ) )
Then, after the last iteration, the messages γ ( h k , n ( j ) ) and ν ( q m , n ( j ) ) are expressed as
γ ( h k , n ( j ) ) = m = 1 M n ( j ) δ m k ( T ) ( h k , n ( j ) )
ν ( q m , n ( j ) ) = k = 1 K n 1 ( j ) ς k m ( T ) ( q m , n ( j ) )
According to γ ( h k , n ( j ) ) , χ ( λ CL , n ( j ) ) and α k , n ( j ) , the message η k ( j ) ( x n ) related to the UE is
η k ( j ) ( x n ) = h k , n ( j ) = 0 M n ( j ) λ CL , n ( j ) G γ ( h k , n ( j ) ) χ ( λ CL , n ( j ) ) b ( x n , a _ k , n ( j ) , 1 , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) × α k ( a _ k , n ( j ) , 1 ) d a _ k , n ( j ) + γ ( h k , n ( j ) = 0 ) α k , n ( j )
5.
Measurement update for legacy PFs
Similar to the measurement update process for the UE, the message about legacy PFs can be expressed as
η ( a _ k , n ( j ) , 1 ) = h k , n ( j ) = 0 M n ( j ) λ CL , n ( j ) G σ n 2 γ ( h k , n ( j ) ) χ ( λ CL , n ( j ) ) b ( x n , a _ k , n ( j ) , 1 , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) α ( x ˜ n , σ n 2 ) d x ˜ n
η k , n ( j ) = η ( a _ k , n ( j ) , 0 ) = γ ( h k , n ( j ) = 0 )
6.
Measurement update for new PFs
The message μ ( a ¯ m , n ( j ) , e ¯ m , n ( j ) ) related to new PFs is calculated as
μ ( a ¯ m , n ( j ) , 1 ) = ν ( q m , n ( j ) = 0 ) χ ( λ CL , n ( j ) ) c ( x n , a ¯ m , n ( j ) , 1 , 0 , λ CL , n ( j ) ; z n ( j ) ) α ( x ˜ n , σ n 2 ) d x ˜ n
μ m , n ( j ) = μ ( a ¯ m , n ( j ) , 0 ) = q m , n ( j ) = 0 K n 1 ( j ) ν ( q m , n ( j ) )
7.
Calculation of beliefs
After calculating the messages above, we can obtain the beliefs for approximating marginal posterior PDFs. First, up to a constant factor, the belief of the UE augmented state f ( x n ) = f ( x ˜ n , σ n 2 ) is expressed as
f ( x ˜ n , σ n 2 ) α ( x ˜ n , σ n 2 ) j = 1 J k K n 1 ( j ) η k ( j ) ( x n )
where the belief f ( x ˜ n , σ n 2 ) provides an approximate estimation of the marginal posterior PDF of the UE augmented state and is used to substitute f ( x ˜ n , σ n 2 | z ) in (26). Furthermore, the belief f _ ( a _ k , n ( j ) , e _ k , n ( j ) ) of the legacy PF augmented states y _ k , n ( j ) [ a _ k , n ( j ) T e _ k , n ( j ) ] T is calculated as
f _ ( a _ k , n ( j ) , 1 ) α k ( a _ k , n ( j ) , 1 ) η ( a _ k , n ( j ) , 1 )
and f _ ( a _ k , n ( j ) , 0 ) α k , n ( j ) η k , n ( j ) . The belief f ¯ ( a ¯ m , n ( j ) , e ¯ m , n ( j ) ) of new PF augmented states y ¯ m , n ( j ) = [ a ¯ m , n ( j ) T e ¯ m , n ( j ) ] T is
f ¯ ( a ¯ m , n ( j ) , 1 ) μ ( a ¯ m , n ( j ) , 1 )
and f ¯ m , n ( j ) = f ¯ ( a ¯ m , n ( j ) , 0 ) μ m , n ( j ) . Then, f _ ( a _ k , n ( j ) , 1 ) and f ¯ ( a ¯ m , n ( j ) , 1 ) are leveraged to approximate the marginal posterior PDF of PF augmented states f ( a k , n ( j ) , e k , n ( j ) = 1 | z ) in (27) jointly, where k K n - 1 ( j ) n ( j ) .
Next, the belief p ( λ CL , n ( j ) ) used to approximate the marginal posterior PMF of the mean number of clutter p ( λ CL , n ( j ) | z ) is expressed as
p ( λ CL , n ( j ) ) = χ ( λ CL , n ( j ) ) β ( λ CL , n ( j ) )
where
β ( λ CL , n ( j ) ) = h k , n ( j ) = 0 M n ( j ) σ n 2 γ ( h k , n ( j ) ) b ( x n , a _ k , n ( j ) , 1 , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) × α k ( a _ k , n ( j ) , 1 ) α ( x ˜ n , σ n 2 ) d a _ k , n ( j ) d x ˜ n + γ ( h k , n ( j ) = 0 ) α k , n ( j )
However, like other nonlinear non-Gaussian filtering models, there are no closed-form solutions for the messages and beliefs above; that is, calculating the beliefs of the UE and PFs directly is infeasible. Therefore, a particle filter-based implementation is used to approximate the beliefs of the UE and PFs (refer to [38,42] for details).

5. Simulation Results and Discussions

This section will evaluate the proposed algorithm from three scenarios and compare it with the conventional BP-SLAM algorithm with unknown parameters and the ideal BP-SLAM algorithm with known parameters.

5.1. Simulation Parameters

We use the environment map shown in Figure 1 as the reference scenario of our simulation. The UE is assumed to follow a nearly constant velocity model, the evolution PDF f ( x n | x n 1 ) of which is x n = A x n 1 + B w n , where A 4 × 4 and B 4 × 2 are defined in [43] (assuming that the sampling period is Δ T = 1 s). The process noise w n is zero-mean and Gaussian with covariance matrix σ p 2 I 2 , where I 2 denotes the identity matrix of size 2 × 2 . Thus, the difference of the UE dynamic models is determined by the UE process noise variance σ p 2 . σ p 2 describes the average increment of the UE speed within the sampling period Δ T , and a σ p 2 with a high value is especially used to model the acceleration and turning process of the UE. Therefore, to accurately describe the motion state of the UE, we propose to use two dynamic models with different process noise variances to represent the motion of the UE. These two process noise variances are expressed as σ p 1 2 and σ p 2 2 , where σ p 1 2 can be chosen as a rule of thumb and σ p 2 2 is at least one order of magnitude larger than σ p 1 2 to describe the state mutation of the UE. Since PFs are static, the evolution model can be considered to be f ( a _ k , n ( j ) | a k , n 1 ( j ) ) = δ ( a _ k , n ( j ) a k , n 1 ( j ) ) , where δ ( ) is the Dirac delta function. However, for the numerical stability of our algorithm, we introduce a small process noise into the evolution model of PFs, i.e., a _ k , n ( j ) = a k , n 1 ( j ) + ϖ k , n ( j ) , where ϖ k , n ( j ) is the zero-mean Gaussian white noise with covariance matrix σ a 2 I 2 .
Different from most SLAM algorithms, our measurement model is only based on range measurements of MPCs in the radio channel, and the range measurement z m , n ( j ) is expressed as
z m , n ( j ) = u n - a k , n ( j ) + e m , n ( j )
where the measurement noise e m , n ( j ) is independent and identically distributed with m, n and j and is a zero-mean Gaussian distribution with variance v m , n ( j ) 2 . We assume that the standard deviation of the measurement noise is v m , n ( j ) =0.1 m. This measurement model (55) determines the factors f ( z m , n ( j ) | x ˜ n , a ¯ m , n ( j ) ) and f ( z m , n ( j ) | x ˜ n , a _ k , n ( j ) ) in the likelihood function (18). In each simulation run, noisy range measurements are produced according to (55). The realization of (55) is based on the UE position u n in its trajectory, the fixed PT positions a 1 ( 1 ) and a 1 ( 2 ) , and the fixed VT positions a l ( j ) 2 , l = 2 , , L n ( j ) , j = 1 , 2 in Figure 1, where L n ( 1 ) = 6 and L n ( 2 ) = 5 are the number of features (including the PT and the associated VTs).
In our proposed algorithm, clutter is modeled as Poisson distribution with mean value λ CL , n ( j ) . We assume that the spatial PDF of clutter is uniform in the region of interest (ROI), which is a circular area with a radius of 30 m around the center of the environment map in Figure 1. Therefore, the clutter intensity can be uniquely determined by the mean value λ CL , n ( j ) . We model λ CL , n ( j ) as a discrete random variable taking values from the finite set G = { 0.1 , 1 , 2 , , 9 } , that is, λ CL , n ( j ) has N g = 10 possible values. The corresponding transition matrix G ( j ) = G is chosen as [33]: when i = 2 , 3 , , 9 , there are [ G ] i , i - 1 = 0.02 , [ G ] i , i = 0.95 and [ G ] i , i + 1 = 0.03 ; furthermore, [ G ] 1 , 1 = 0.99 , [ G ] 1 , 2 = 0.01 , [ G ] 10 , 9 = 0.1 and [ G ] 10 , 10 = 0.9 ; otherwise, [ G ] i , j = 0 . The PDF f new , n ( a ¯ m , n ( j ) ) and mean value λ new , n ( j ) of the newly detected features are estimated online by the PHD filter, which uses the spatial PDF f b , n ( a ¯ m , n ( j ) ) and mean number λ b ( j ) of newborn features [44]. We assume that f b , n ( a ¯ m , n ( j ) ) is uniform in the ROI. When n = 1, the initialization of newly detected features is realized by setting λ new , 1 ( j ) = 6 and the initial PDF f new , 1 ( a ¯ m , 1 ( j ) ) to be uniform in the ROI. In addition, the detection probability P d ( j ) ( x n , y k , n ( j ) ) is considered fixed, i.e., P d ( j ) ( x n , y k , n ( j ) ) = P d .
The messages and beliefs in Section 4.3 are approximately calculated via the particle filter-based implementation with N par = 10000 particles in the proposed algorithm. Particles representing the initial state of the PT follow the two-dimensional Gaussian distribution N ( a 1 , 1 ( j ) , σ b , 1 2 I 2 ) , where a 1 , 1 ( j ) denotes the initial positions of PT1 (j = 1) and PT2 (j = 2). Particles representing the initial state of the UE follow a 4-dimensional uniform distribution with the center x 1 = [ p 1 T 0 0 ] T , where p 1 is the starting position of the UE. This implies that the prior information of the initial state of the PT and the UE is known, but the prior information related to VTs is unknown. Finally, the termination condition of the iterative DA is [ k = 1 k = K n 1 ( j ) m = 1 m = M n ( j ) ( δ m k ( i ) ( h k , n ( j ) ) δ m k ( i 1 ) ( h k , n ( j ) ) ) 2 ] 1 / 2 < 10 7 or the maximum number of iterations I max = 1000 is reached. Next, we will divide the scenarios into groups to evaluate the performance of our proposed algorithm and provide the corresponding numerical results. The common simulation parameters of all scenarios are shown in Table 1. The parameter settings in Table 1 are based on [25] and it is verified in [25] that the BP-SLAM algorithm can achieve excellent performance under these conditions.

5.2. First Scenario: Unknown Process Noise Only

In the first scenario, the UE switches between two nearly constant velocity models with different process noise variances, σ p 1 2 = 0.01 2 and σ p 2 2 = 0.1 2 . The simulation time is 800 s. For 0–720 s, the UE process noise variance is σ p 1 2 = 0.01 2 ; then, the UE accelerates suddenly with the process noise variance changing to σ p 2 2 = 0.1 2 for 720–730 s; finally, for 730–800 s, the UE process noise variance is the same as that for the first stage. The transition matrix of the UE process noise variance is [ Q b ] 1 , 1 = [ Q b ] 2 , 2 = 0.9975 and [ Q b ] 1 , 2 = [ Q b ] 2 , 1 = 0.0025 . In this scenario, the clutter intensity is considered to be a fixed prior. We assume that the mean number of clutter is λ CL , n ( j ) = 0 .1 . We compare our proposed algorithm with the conventional BP-SLAM algorithm [25] with unknown process noise variance and the ideal BP-SLAM algorithm with known process noise variance, which can obtain the prior information of the UE process noise variance at each time in advance. The ideal BP-SLAM algorithm is impossible in practice but can be used as a benchmark to evaluate the performance of our proposed algorithm. In the conventional BP-SLAM algorithm with unknown process noise variance, we considered two different parameter settings labeled BP-SLAM A and BP-SLAM B. The UE process noise variances in BP-SLAM A and BP-SLAM B both remain unchanged, but their settings are different, i.e., σ p 2 = σ p 1 2 = 0.01 2 for BP-SLAM A and σ p 2 = σ p 2 2 = 0.1 2 for BP-SLAM B.
Figure 3 shows the performance of the proposed algorithm in the first scenario, which is obtained based on the average of 100 simulation runs. Figure 3a,b depict the mean optimal subpattern assignment (MOSPA) [45] errors for PT1 and PT2 and their associated VTs, respectively. The MOSPA error is based on Euclidean distance with a cutoff parameter of 5 m and order 1. The MOSPA error combines the estimation errors of correctly detected PFs and incorrectly detected PFs, which is often used as the standard to evaluate the mapping accuracy of multipath-assisted SLAM algorithms. In Figure 3a, period n [ 0 , 720 ] , the MOSPA error of BP-SLAM B for PT1 and its associated VTs is always higher than the other three algorithms. Then, when the UE process noise variance suddenly changes, i.e., n [ 720 , 730 ] , the MOSPA error of BP-SLAM A rises sharply with a maximum of more than 3 m. However, the MOSPA error of the proposed algorithm is always smaller than that of BP-SLAM A and BP-SLAM B in this period. This relationship still exists even after the UE returns to the normal motion state, i.e., n [ 730 , 800 ] . The MOSPA error for PT2 and its associated VTs shown in Figure 3b also have similar performance. Figure 3c,d describe the average number of detected PFs associated with PT1 and PT2, respectively. From these two figures, we can see that the average number of detected PFs for the proposed algorithm is always close to the actual number of features in the environment, namely, L n ( 1 ) = 6 and L n ( 2 ) = 5 . However, the average number of detected PFs for BP-SLAM A deviates from the truth after 720 s, i.e., the UE process noise variance changes suddenly. The average number of detected PFs for BP-SLAM B deviates from the truth even when the UE process noise remains unchanged, i.e., n [ 0 , 720 ] . Figure 3e depicts the root mean square errors (RMSEs) of the UE position for the four algorithms. Figure 3e shows that in period n [ 720 , 730 ] , the UE position RMSE for BP-SLAM A rise sharply with a maximum of more than 1.2 m, while the UE position RMSE of the proposed algorithm is always less than 0.4 m in this period, which is close to the ideal BP-SLAM algorithm with known parameters. Although the UE position RMSE for BP-SLAM B has no sharp rise in this period, it is higher than that of the proposed algorithm all the time.
Table 2 summarizes the error statistics of different algorithms in the first scenario within 800 s, including the average UE position RMSE, the average MOSPA error and the average error for the estimated number of detected PFs. The latter two errors are based on the average of PT1 and PT2 statistics. From Table 2, we can see that the statistical error of the proposed algorithm is significantly smaller than that of BP-SLAM A and BP-SLAM B, and close to the ideal BP-SLAM algorithm with known parameters. Therefore, Figure 3 and Table 2 both illustrate that the proposed algorithm is robust to the uncertainty of the process noise.

5.3. Second Scenario: Unknown Clutter Intensity Only

In the second scenario, we assume that only the clutter intensity parameter is unknown, which means that the UE does not have a sudden change in the dynamic model, with the process noise variance keeping σ p 2 = 0.01 2 during the movement. Therefore, it is unnecessary to consider the uncertainty of process noise in this scenario. The simulation time is also 800 s. Clutter is generated by a Poisson distribution with mean number λ CL , n ( j ) = 6 , and the spatial PDF of clutter is uniform in the ROI, i.e., f CL ( z m , n ( j ) ) = 1 / D , where D is the area of the ROI. Thus, the clutter intensity can be expressed as 6 / D . We also compare the proposed algorithm with the conventional BP-SLAM algorithm with unknown parameters and the ideal BP-SLAM algorithm with known parameters in this scenario. Similar to the first scenario, we consider two parameter settings in the conventional BP-SLAM algorithm labeled BP-SLAM A and BP-SLAM B. They have no prior information on the true clutter intensity and just set this parameter by experience. BP-SLAM A and BP-SLAM B initialize the clutter intensity to 0.1 / D and 12 / D , respectively. The ideal BP-SLAM algorithm with known parameters can obtain the prior information of the true clutter intensity in advance, which can also be used as the benchmark to evaluate the performance of our proposed algorithm in this scenario.
Figure 4 illustrates the simulation results of the proposed algorithm in the second scenario, which are also obtained based on the average of 100 simulation runs. Figure 4a,b show the MOSPA errors for PT1 and PT2 and their associated VTs, respectively. From Figure 4a, we see that the MOSPA errors of the four algorithms for PT1 all have a decreasing trend and finally reach steady states. The MOSPA errors for PT1 and the associated VTs of BP-SLAM A and BP-SLAM B converge at approximately 0.6 m and 0.4 m, respectively, while the performance of the proposed algorithm is better, converging near 0.2 m, which is close to the result of the BP-SLAM algorithm with known parameters. The results in Figure 4b are similar to those of Figure 4a. When the MOSPA error of the four algorithms becomes steady, the MOSPA error of the proposed algorithm is significantly smaller than that of BP-SLAM A and BP-SLAM B. Figure 4c,d describe the average number of detected PFs associated with PT1 and PT2 every 60 s, respectively. These two figures show that the average numbers of detected PFs by the proposed algorithm converge to the true values within 60 s, while the numbers of detected PFs by BP-SLAM A are always in a disordered state. Meanwhile, although the results for BP-SLAM B are convergent, they converge to the wrong number. Figure 4e depicts the UE position RMSEs for the four algorithms. Figure 4e shows that the position accuracy of the proposed algorithm is significantly better than that of BP-SLAM A and BP-SLAM B. During the simulation run time of 800 s, the RMSE of the proposed algorithm is always lower than 0.15 m, which is close to the performance of the BP-SLAM algorithm with known parameters. However, the UE position RMSE of BP-SLAM A increases rapidly in 800 s, and the maximum is more than 0.65 m. Although the performance of BP-SLAM B is relatively good, its position accuracy is lower than that of the proposed algorithm all the time.
Table 3 summarizes the error statistics of different algorithms in the second scenario within 800 s. The MOSPA error and the error for the estimated number of detected PFs are also based on the average of PT1 and PT2 statistics. It can be seen from Table 3 that the statistical error of the proposed algorithm is significantly smaller than that of BP-SLAM A and BP-SLAM B, and close to the ideal BP-SLAM algorithm with known parameters. Therefore, Figure 4 and Table 3 both verify that the proposed algorithm is robust to the uncertainty of the clutter intensity.

5.4. Third Scenario: Unknown Process Noise and Clutter Intensity

In the third scenario, the uncertainty of the process noise and the clutter intensity are considered together, and the simulation time is 800 s. The variation in the UE process noise variance is the same as that in the first scenario. Most of the time, the UE process noise variance is σ p 1 2 = 0.01 2 , but when n [ 720 , 730 ] , due to the sudden acceleration, the UE process noise variance changes to σ p 2 2 = 0.1 2 . The BP-SLAM algorithm with unknown parameters considers that the UE process noise variance is fixed to σ p 2 = σ p 1 2 = 0.01 2 , while our proposed algorithm can adaptively switch between the two dynamic models with the transition matrix of the process noise variance the same as that in the first scenario. The generation mechanism of clutter is consistent with the second scenario. The true clutter intensity is 4 / D , but the proposed algorithm and the BP-SLAM algorithm with unknown parameters do not have this prior information. We also consider two parameter settings in the conventional BP-SLAM algorithm with unknown parameters in this scenario, labeled BP-SLAM A and BP-SLAM B. In BP-SLAM A, the UE process noise variance is fixed to σ p 2 = σ p 1 2 = 0.01 2 and the clutter intensity is initialized to 0.1 / D by experience. In BP-SLAM B, these two parameters are set to σ p 2 = σ p 2 2 = 0.1 2 and 8 / D , respectively. In the third scenario, the BP-SLAM algorithm with known parameters can obtain the prior parameters of the process noise and the clutter intensity at each time in advance, which can still be used as the benchmark to evaluate the performance of the proposed algorithm.
Figure 5 depicts the performance of the proposed algorithm in the third scenario, which is also obtained based on the average of 100 simulation runs. Figure 5a,b describe the MOSPA errors for PT1 and PT2 and their associated VTs, respectively. Figure 5a shows that the MOSPA errors of the four algorithms all reach steady states within approximately 150 s, and BP-SLAM A and BP-SLAM B are stable at approximately 1.8 m and 0.5 m, respectively. However, the performance of the proposed algorithm and the BP-SLAM algorithm with known parameters is obviously better with both below 0.3 m. In particular, the MOSPA error of BP-SLAM A rises sharply during the sudden change of the UE process noise variance, i.e., n [ 720 , 730 ] , and the maximum arrives near 3.2 m. The MOSPA errors of the proposed algorithm and the BP-SLAM algorithm with known parameters do not change significantly in this period. The MOSPA error for PT2 and the associated VTs described in Figure 5b has a similar performance to that in Figure 5a, i.e., regardless of whether the UE process noise changes, the performance of the proposed algorithm is better than that of BP-SLAM A and BP-SLAM B. Figure 5c,d show the average number of detected PFs associated with PT1 and PT2 every 60 s, respectively. These two figures show that the average numbers of detected PFs by the proposed algorithm are always near the true value, while the results of BP-SLAM A are not convergent. Although the results for BP-SLAM B are convergent, they converge on the wrong number. Therefore, in the third scenario, the estimation performance of the proposed algorithm for the average number of detected PFs is obviously better than that of BP-SLAM A and BP-SLAM B. Figure 5e describes the UE position RMSEs for the four algorithms in the third scenario. The time average RMSE of the proposed algorithm is 0.11 m within 800 s, while that of BP-SLAM A and BP-SLAM B are 1.46 m and 0.25 m, respectively. In particular, the time average RMSE of the proposed algorithm is 0.37 m when the UE process noise variance changes suddenly, i.e., n [ 720 , 730 ] , while BP-SLAM A and BP-SLAM B are 3.18 m and 0.55 m in this period. For the BP-SLAM algorithm with known parameters, the time average RMSE of the UE position is 0.07 m during 800 s and 0.11 m when n [ 720 , 730 ] . Based on the above data, the UE position accuracy of the proposed algorithm is significantly better than that of the BP-SLAM algorithm with unknown parameters and close to that of the BP-SLAM algorithm with known parameters.
Table 4 summarizes the error statistics of different algorithms in the third scenario within 800 s. The MOSPA error and the error for the estimated number of detected PFs are also based on the average of PT1 and PT2 statistics. From Table 4, we can see that the statistical error of the proposed algorithm is significantly smaller than that of BP-SLAM A and BP-SLAM B, and close to the ideal BP-SLAM algorithm with known parameters. Therefore, Figure 5 and Table 4 show that the proposed algorithm is robust to both the clutter intensity and the process noise uncertainty.

5.5. Calculation Complexity Analysis

The computational complexity of the conventional BP-SLAM algorithm is proportional to O ( N par ) , where N par is the number of particles [25]. The proposed algorithm adds the evolutions of the process noise variance and the clutter intensity to the Bayesian model compared with the conventional BP-SLAM algorithm; hence, its computational complexity is proportional to O ( N par + N b 2 + N g 2 ) [33], where N b and N g denote the number of process noise variances and clutter intensity, respectively. Considering that N b and N g are generally small compared to the number of particles N par , and therefore our proposed algorithm has nearly the same computational complexity as the conventional BP-SLAM algorithm. The average CPU run time consumed by each time step n in the third scenario is also used to compare the calculation complexity of different algorithms, where N par = 10000 , N b = 2 and N g = 10 . The results of CPU run-time are based on 100 simulation runs and 800 time steps using a MATLAB implementation on an Intel Core i7-3632QM [email protected] GHz. Table 5 shows the complexity comparison of different algorithms in terms of both the complexity orders and the CPU run times. From Table 5, we can see that the proposed algorithm has almost the same calculation complexity as the conventional BP-SLAM algorithm with unknown parameters, but the performance of positioning and mapping accuracy is much better, close to the ideal BP-SLAM algorithm with known parameters. The reason is that the proposed algorithm considers the influence of the process noise and the clutter intensity uncertainty.

6. Conclusions

In multipath-assisted SLAM using radio signals, the estimation of the process noise and the clutter intensity is difficult in practice. Therefore, it is highly significant to study the multipath-assisted SLAM algorithm with the ability to adaptively estimate these unknown parameters. A robust multipath-assisted SLAM algorithm is proposed in this paper to solve the problem that the accuracy of the multipath-assisted SLAM algorithm decreases when the process noise and the clutter intensity are statistically uncertain. The proposed algorithm can still work with unknown time-varying parameters by continuously inferring an unknown parameter model, the states of the UE and environmental features. Specifically, we integrate the process noise variance and the clutter intensity into the factor graph representing the statistical structure of the SLAM problem and use the Markov chain model to describe the evolution of these unknown parameters. Finally, the BP message passing algorithm is performed to calculate the marginal posterior PDFs of each parameter to achieve the goal of simultaneous localization and mapping and adaptively learning the process noise variance and the clutter intensity.
The simulation results demonstrate that the proposed algorithm is superior to the conventional BP-SLAM algorithm in positioning and mapping accuracy with the uncertainty of the process noise and the clutter intensity. In particular, when the process noise and the clutter intensity are both unknown, the UE position time average RMSE at all times of the proposed algorithm is significantly less than that of the conventional BP-SLAM algorithm (0.11 m versus 1.46 m), and the performance in the MOSPA error and the average number of detected PFs is also obviously better than that of the conventional BP-SLAM algorithm (0.77 m versus 2.24 m and 0.06 versus 8.5, respectively), taking BP-SLAM A as an example. Therefore, the proposed algorithm can effectively restrain the impact of the process noise and the clutter intensity uncertainty on SLAM estimation accuracy and improve the robustness of a multipath-assisted SLAM system using radio signals.
We assume time synchronization between the transmitters and the UE in this paper. A promising direction for future research is the extension of our robust multipath-assisted SLAM algorithm to an unsynchronized sensor network.

Author Contributions

Methodology, Z.D.; Project administration, B.L.; Software, Z.D.; Writing—Original draft, Z.D.; Writing—Review & Editing, B.L. and C.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Natural Science Basic Research Program of Shaanxi, grant number 2021JQ-122.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Durrant-Whyte, H.; Bailey, T. Simultaneous localization and mapping: Part I. IEEE Robot. Autom. Mag. 2006, 13, 99–110. [Google Scholar] [CrossRef] [Green Version]
  2. Bailey, T.; Durrant-Whyte, H. Simultaneous localization and mapping (SLAM): Part II. IEEE Robot. Autom. Mag. 2006, 13, 108–117. [Google Scholar] [CrossRef] [Green Version]
  3. Lajoie, P.Y.; Ramtoula, B.; Chang, Y. DOOR-SLAM: Distributed, Online, and Outlier Resilient SLAM for Robotic Teams. IEEE Robot. Autom. Lett. 2020, 5, 1656–1663. [Google Scholar] [CrossRef] [Green Version]
  4. Chiang, K.W.; Tsai, G.J.; Chu, H.J. Performance Enhancement of INS/GNSS/Refreshed-SLAM Integration for Acceptable Lane-Level Navigation Accuracy. IEEE Trans. Veh. Technol. 2020, 69, 2463–2476. [Google Scholar] [CrossRef]
  5. Mendrzik, R.; Meyer, F.; Bauch, G. Enabling Situational Awareness in Millimeter Wave Massive MIMO Systems. IEEE J. Sel. Top. Signal Process. 2019, 13, 1196–1211. [Google Scholar] [CrossRef]
  6. Chen, Z.; Zhu, G.; Wang, S. M3: Multipath Assisted Wi-Fi Localization with a Single Access Point. IEEE Trans. Mob. Comput. 2021, 20, 588–602. [Google Scholar] [CrossRef]
  7. Ulmschneider, M.; Gentner, C.; Jost, T. Association of Transmitters in Multipath-Assisted Positioning. In Proceedings of the IEEE Global Communications Conference (GLOBECOM), Singapore, 4–8 December 2017. [Google Scholar]
  8. Li, Z.; Tian, Z.; Wang, Z. Multipath-Assisted Indoor Localization Using a Single Receiver. IEEE Sens. J. 2021, 21, 692–705. [Google Scholar] [CrossRef]
  9. Li, Z.; Tian, Z.; Wang, Z. Multipath-Assisted Indoor Localization: Turning Multipath Signal from Enemy to Friend. In Proceedings of the IEEE Global Communications Conference (GLOBECOM), Waikoloa, HI, USA, 9–13 December 2019. [Google Scholar]
  10. Wielandt, S.; Strycker, L.D. Indoor Multipath Assisted Angle of Arrival Localization. Sensors 2017, 17, 2522. [Google Scholar] [CrossRef] [Green Version]
  11. Witrisal, K.; Meissner, P.; Leitinger, E. High-accuracy localization for assisted living: 5G systems will turn multipath channels from foe to friend. IEEE Signal Process. Mag. 2016, 33, 59–70. [Google Scholar] [CrossRef]
  12. Guidi, F.; Closas, P.; Dardari, D. Personal mobile radars with millimeter-wave massive arrays for indoor mapping. IEEE Trans. Mobile Comput. 2016, 15, 1471–1484. [Google Scholar] [CrossRef]
  13. Guerra, A.; Guidi, F.; Dardari, D. Occupancy grid mapping for personal radar applications. In Proceedings of the IEEE Statistical Signal Processing Workshop (SSP), Freiburg, Germany, 10–13 June 2018. [Google Scholar]
  14. Zhang, H.; Tan, S.Y. TOA based indoor localization and tracking via single-cluster PHD filtering. In Proceedings of the IEEE Global Communications Conference (GLOBECOM), Singapore, 4–8 December 2017. [Google Scholar]
  15. Kim, H.; Granström, K.; Gao, L. 5G mmWave Cooperative Positioning and Mapping Using Multi-Model PHD Filter and Map Fusion. IEEE Trans. Wirel. Commun. 2020, 19, 3782–3795. [Google Scholar] [CrossRef] [Green Version]
  16. Mendrzik, R.; Wymeersch, H.; Bauch, G. Joint Localization and Mapping through Millimeter Wave MIMO in 5G Systems. In Proceedings of the IEEE Global Communications Conference (GLOBECOM), Abu Dhabi, United Arab Emirates, 9–13 December 2018. [Google Scholar]
  17. Leitinger, E.; Meissner, P.; Rudisser, C. Evaluation of position-related information in multipath components for indoor positioning. IEEE J. Sel. Areas Commun. 2015, 33, 2313–2328. [Google Scholar] [CrossRef] [Green Version]
  18. Leitinger, E.; Meissner, P.; Witrisal, K. Simultaneous localization and mapping using multipath channel information. In Proceedings of the IEEE International Conference on Communication Workshop (ICCW), London, UK, 8–12 June 2015. [Google Scholar]
  19. Gentner, C.; Ma, B.; Ulmschneider, M. Simultaneous localization and mapping in multipath environments. In Proceedings of the IEEE/ION Position, Location and Navigation Symposium (PLANS), Savannah, GA, USA, 11–14 April 2016. [Google Scholar]
  20. Li, Z.; Tian, Z.; Wang, Z. WalkAround: Multipath-assisted Indoor Localization and Mapping Using a Single Receiver. In Proceedings of the IEEE International Conference on Communications (ICC), Dublin, Ireland, 7–11 June 2020. [Google Scholar]
  21. Ulmschneider, M.; Gentner, C. Multipath assisted positioning for pedestrians using LTE signals. In Proceedings of the IEEE/ION Position, Location and Navigation Symposium (PLANS), Savannah, GA, USA, 11–14 April 2016. [Google Scholar]
  22. Gentner, C.; Jost, T.; Wang, W. Multipath Assisted Positioning with Simultaneous Localization and Mapping. IEEE Trans. Wirel. Commun. 2016, 15, 6104–6117. [Google Scholar] [CrossRef] [Green Version]
  23. Karásek, R.; Gentner, C. Stochastic Data Association for Multipath Assisted Positioning Using a Single Transmitter. IEEE Access 2020, 8, 46735–46752. [Google Scholar] [CrossRef]
  24. Liu, Y.; Lian, B.; Zhou, T. Gaussian message pasing-based cooperative localization with node selection scheme in wireless networks. Signal Process. 2019, 156, 166–176. [Google Scholar] [CrossRef]
  25. Leitinger, E.; Meyer, F.; Hlawatsch, F. A Belief Propagation Algorithm for Multipath-Based SLAM. IEEE Trans. Wirel. Commun. 2019, 18, 5613–5629. [Google Scholar] [CrossRef] [Green Version]
  26. Leitinger, E.; Meyer, F.; Tufvesson, F. Factor Graph Based Simultaneous Localization and Mapping Using Multipath Channel Information. In Proceedings of the IEEE International Conference on Communications Workshops (ICC Workshops), Paris, France, 21–25 May 2017. [Google Scholar]
  27. Leitinger, E.; Grebien, S.; Li, X. On the Use of Mpc Amplitude Information in Radio Signal Based SLAM. In Proceedings of the IEEE Statistical Signal Processing Workshop (SSP), Freiburg, Germany, 10–13 June 2018. [Google Scholar]
  28. Leitinger, E.; Grebien, S.; Witrisal, K. Multipath-based SLAM Exploiting AoA and Amplitude Information. In Proceedings of the IEEE International Conference on Communications Workshops (ICC Workshops), Shanghai, China, 20–24 May 2019. [Google Scholar]
  29. Bar-Shalom, Y.; Li, X.-R. Multitarget-Multisensor Tracking: Principles and Techniques; Yaakov Bar-Shalom: Storrs, CT, USA, 1995; pp. 1–4. [Google Scholar]
  30. Mahler, R.P.S.; Vo, B.T.; Vo, B.N. CPHD Filtering with Unknown Clutter Rate and Detection Profile. IEEE Trans. Signal. Process. 2011, 59, 3497–3513. [Google Scholar] [CrossRef]
  31. Vo, B.T.; Vo, B.N.; Hoseinnezhad, R. Robust multi-Bernoulli filtering. IEEE J. Sel. Top. Signal Process. 2013, 7, 399–409. [Google Scholar] [CrossRef]
  32. Ristic, B.; Clark, D.; Vo, B.T.; Vo, B.N. Adaptive target birth intensity for PHD and CPHD filters. IEEE Trans. Aerosp. Electron. Syst. 2012, 48, 1656–1668. [Google Scholar] [CrossRef]
  33. Soldi, G.; Meyer, F.; Braca, P.; Hlawatsch, F. Self-Tuning Algorithms for Multisensor-Multitarget Tracking Using Belief Propagation. IEEE Trans. Signal. Process. 2019, 67, 3922–3937. [Google Scholar] [CrossRef]
  34. Kschischang, F.R.; Frey, B.J.; Loeliger, H.-A. Factor graphs and the sum-product algorithm. IEEE Trans. Inf. Theory 2001, 47, 498–519. [Google Scholar] [CrossRef] [Green Version]
  35. Salmi, J.; Richter, A.; Koivunen, V. Detection and tracking of MIMO propagation path parameters using state-space approach. IEEE Trans. Signal. Process. 2009, 57, 1538–1550. [Google Scholar] [CrossRef] [Green Version]
  36. Mahler, R.P.S. Statistical Multisource-Multitarget Information Fusion; Artech House: Norwood, MA, USA, 2007. [Google Scholar]
  37. Jo, K.; Chu, K.; Sunwoo, M. Interacting Multiple Model Filter-Based Sensor Fusion of GPS With In-Vehicle Sensors for Real-Time Vehicle Positioning. IEEE Trans. Intell. Transp. 2012, 13, 329–343. [Google Scholar] [CrossRef]
  38. Meyer, F.; Braca, P.; Willett, P.; Hlawatsch, F. A scalable algorithm for tracking an unknown number of targets using multiple sensors. IEEE Trans. Signal. Process. 2017, 65, 3478–3493. [Google Scholar] [CrossRef] [Green Version]
  39. Wymeersch, H.; Seco-Granados, G. Adaptive Detection Probability for mmWave 5G SLAM. In Proceedings of the IEEE 2nd 6G Wireless Summit (6G SUMMIT), Levi, Finland, 17–20 March 2020. [Google Scholar]
  40. Williams, J.L.; Lau, R. Approximate evaluation of marginal association probabilities with belief propagation. IEEE Trans. Aerosp. Electron. Syst. 2014, 50, 2942–2959. [Google Scholar] [CrossRef] [Green Version]
  41. Meyer, F.; Kropfreiter, T.; Williams, J.L. Message passing algorithms for scalable multitarget tracking. Proc. IEEE 2018, 106, 221–259. [Google Scholar] [CrossRef]
  42. Meyer, F.; Hlinka, O.; Wymeersch, H.; Hlawatsch, F. Distributed localization and tracking of mobile networks including noncooperative objects. IEEE Trans. Signal Inf. Process. Netw. 2016, 2, 57–71. [Google Scholar] [CrossRef] [Green Version]
  43. Li, X.R.; Jilkov, V.P. Survey of maneuvering target tracking. part I. dynamic models. IEEE Trans. Aerosp. Electron. Syst. 2003, 39, 1333–1364. [Google Scholar] [CrossRef]
  44. Horridge, P.; Maskell, S. Using a probabilistic hypothesis density filter to confirm tracks in a multi-target environment. In Proceedings of the Informatik Schafft Communities (INFORMATIK-11), Berlin, Germany, 4–7 July 2011. [Google Scholar]
  45. Schuhmacher, D.; Vo, B.T.; Vo, B.N. A consistent metric for performance evaluation of multi-object filters. IEEE Trans. Signal Process. 2008, 56, 3447–3457. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Example of an environment map.
Figure 1. Example of an environment map.
Remotesensing 13 01625 g001
Figure 2. Factor graph representing the factorization of the joint posterior probability density function (PDF) in (30). All factor nodes, variable nodes and messages related to the user equipment (UE) augmented state are represented in red color, those related to legacy potential features (PF) states are represented in purple color, those related to new PF states are represented in blue color, and the section of the clutter intensity and the iterative data association (DA) are represented in orange color and green color, respectively. The abbreviations are defined as K K n ( j ) , M M n ( j ) , x x n , y _ k y _ k , n ( j ) , h k h k , n ( j ) , q m q m , n ( j ) , λ CL λ CL , n ( j ) , φ k , m φ ( h k , n ( j ) , q m , n ( j ) ) , c m c ( x n , a ¯ m , n ( j ) , e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) , b k b ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) , f f ( x n | x n 1 ) , f k f ( y _ k , n ( j ) | y k , n 1 ( j ) ) , p p ( λ CL , n ( j ) | λ CL , n 1 ( j ) ) , α α ( x n ) , α k α k ( a _ k , n ( j ) , e _ k , n ( j ) ) , f f ( x n 1 ) , f f ( x n ) , f _ k f _ ( a k , n 1 ( j ) , e k , n 1 ( j ) ) , f _ k f _ ( a _ k , n ( j ) , e _ k , n ( j ) ) , f ¯ m f ¯ ( a ¯ m , n ( j ) , e ¯ m , n ( j ) ) , p p ( λ CL , n 1 ( j ) ) , p p ( λ CL , n ( j ) ) , β k β ( λ CL , n ( j ) ) , χ χ ( λ CL , n ( j ) ) , ε k ε ( h k , n ( j ) ) , ξ m ξ ( q m , n ( j ) ) , γ k γ ( h k , n ( j ) ) , ν m ν ( q m , n ( j ) ) , δ m , k δ m k ( t ) ( h k , n ( j ) ) , ς k , m ς k m ( t ) ( q m , n ( j ) ) , η k ( j ) η k ( j ) ( x n ) , η k η ( a _ k , n ( j ) , e _ k , n ( j ) ) , and μ m μ ( a ¯ m , n ( j ) , e ¯ m , n ( j ) ) .
Figure 2. Factor graph representing the factorization of the joint posterior probability density function (PDF) in (30). All factor nodes, variable nodes and messages related to the user equipment (UE) augmented state are represented in red color, those related to legacy potential features (PF) states are represented in purple color, those related to new PF states are represented in blue color, and the section of the clutter intensity and the iterative data association (DA) are represented in orange color and green color, respectively. The abbreviations are defined as K K n ( j ) , M M n ( j ) , x x n , y _ k y _ k , n ( j ) , h k h k , n ( j ) , q m q m , n ( j ) , λ CL λ CL , n ( j ) , φ k , m φ ( h k , n ( j ) , q m , n ( j ) ) , c m c ( x n , a ¯ m , n ( j ) , e ¯ m , n ( j ) , q m , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) , b k b ( x n , a _ k , n ( j ) , e _ k , n ( j ) , h k , n ( j ) , λ CL , n ( j ) ; z n ( j ) ) , f f ( x n | x n 1 ) , f k f ( y _ k , n ( j ) | y k , n 1 ( j ) ) , p p ( λ CL , n ( j ) | λ CL , n 1 ( j ) ) , α α ( x n ) , α k α k ( a _ k , n ( j ) , e _ k , n ( j ) ) , f f ( x n 1 ) , f f ( x n ) , f _ k f _ ( a k , n 1 ( j ) , e k , n 1 ( j ) ) , f _ k f _ ( a _ k , n ( j ) , e _ k , n ( j ) ) , f ¯ m f ¯ ( a ¯ m , n ( j ) , e ¯ m , n ( j ) ) , p p ( λ CL , n 1 ( j ) ) , p p ( λ CL , n ( j ) ) , β k β ( λ CL , n ( j ) ) , χ χ ( λ CL , n ( j ) ) , ε k ε ( h k , n ( j ) ) , ξ m ξ ( q m , n ( j ) ) , γ k γ ( h k , n ( j ) ) , ν m ν ( q m , n ( j ) ) , δ m , k δ m k ( t ) ( h k , n ( j ) ) , ς k , m ς k m ( t ) ( q m , n ( j ) ) , η k ( j ) η k ( j ) ( x n ) , η k η ( a _ k , n ( j ) , e _ k , n ( j ) ) , and μ m μ ( a ¯ m , n ( j ) , e ¯ m , n ( j ) ) .
Remotesensing 13 01625 g002
Figure 3. Results for the first scenario: (a) mean optimal subpattern assignment (MOSPA) error for physical transmitter (PT) 1 and the associated virtual transmitters (VTs); (b) MOSPA error for PT2 and the associated VTs; (c) average number of detected PFs associated with PT1 every 60 s; (d) average number of detected PFs associated with PT2 every 60 s; and (e) user equipment (UE) position root mean square errors (RMSEs).
Figure 3. Results for the first scenario: (a) mean optimal subpattern assignment (MOSPA) error for physical transmitter (PT) 1 and the associated virtual transmitters (VTs); (b) MOSPA error for PT2 and the associated VTs; (c) average number of detected PFs associated with PT1 every 60 s; (d) average number of detected PFs associated with PT2 every 60 s; and (e) user equipment (UE) position root mean square errors (RMSEs).
Remotesensing 13 01625 g003
Figure 4. Results for the second scenario: (a) MOSPA error for PT1 and the associated VTs; (b) MOSPA error for PT2 and the associated VTs; (c) average number of detected PFs associated with PT1 every 60 s; (d) average number of detected PFs associated with PT2 every 60 s; and (e) UE position RMSE.
Figure 4. Results for the second scenario: (a) MOSPA error for PT1 and the associated VTs; (b) MOSPA error for PT2 and the associated VTs; (c) average number of detected PFs associated with PT1 every 60 s; (d) average number of detected PFs associated with PT2 every 60 s; and (e) UE position RMSE.
Remotesensing 13 01625 g004
Figure 5. Results for the third scenario: (a) MOSPA error for PT1 and the associated VTs; (b) MOSPA error for PT2 and the associated VTs; (c) average number of detected PFs associated with PT1 every 60 s; (d) average number of detected PFs associated with PT2 every 60 s; and (e) UE position RMSE.
Figure 5. Results for the third scenario: (a) MOSPA error for PT1 and the associated VTs; (b) MOSPA error for PT2 and the associated VTs; (c) average number of detected PFs associated with PT1 every 60 s; (d) average number of detected PFs associated with PT2 every 60 s; and (e) UE position RMSE.
Remotesensing 13 01625 g005
Table 1. Parameter setup.
Table 1. Parameter setup.
σ a σ b , 1 λ n e w , 1 ( j ) v m , n ( j ) P d λ b ( j ) P s P t h P p r N p a r Simulation Runs
10−4 m10−3 m60.1 m0.9510−40.9990.510−410,000100
Table 2. Error statistics of different algorithms in the first scenario.
Table 2. Error statistics of different algorithms in the first scenario.
AlgorithmUE RMSEMOSPANo. of PFs
BP-SLAM A0.2891 m1.1143 m0.1090
BP-SLAM B0.1984 m1.0442 m0.2564
Proposed algorithm0.1605 m0.8759 m0.0005
BP-SLAM (known params)0.0949 m0.7739 m0 (True)
Table 3. Error statistics of different algorithms in the second scenario.
Table 3. Error statistics of different algorithms in the second scenario.
AlgorithmUE RMSEMOSPANo. of PFs
BP-SLAM A0.3768 m1.2199 m8.0151
BP-SLAM B0.1083 m1.0006 m0.8329
Proposed algorithm0.0741 m0.7659 m0.0678
BP-SLAM (known params)0.0714 m0.7376 m0 (True)
Table 4. Error statistics of different algorithms in the third scenario.
Table 4. Error statistics of different algorithms in the third scenario.
AlgorithmUE RMSEMOSPANo. of PFs
BP-SLAM A1.4636 m2.2396 m8.4696
BP-SLAM B0.2452 m1.0713 m0.8099
Proposed algorithm0.1121 m0.7659 m0.0628
BP-SLAM (known params)0.0699 m0.6738 m0 (True)
Table 5. Complexity comparison of different algorithms.
Table 5. Complexity comparison of different algorithms.
AlgorithmComplexity OrderCPU Run Time
BP-SLAM (unknown params) O ( N par ) 0.0645 s
Proposed algorithm O ( N par + N b 2 + N g 2 ) 0.0761 s
BP-SLAM (known params) O ( N par ) 0.0652 s
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Dan, Z.; Lian, B.; Tang, C. Robust Multipath-Assisted SLAM with Unknown Process Noise and Clutter Intensity. Remote Sens. 2021, 13, 1625. https://doi.org/10.3390/rs13091625

AMA Style

Dan Z, Lian B, Tang C. Robust Multipath-Assisted SLAM with Unknown Process Noise and Clutter Intensity. Remote Sensing. 2021; 13(9):1625. https://doi.org/10.3390/rs13091625

Chicago/Turabian Style

Dan, Zesheng, Baowang Lian, and Chengkai Tang. 2021. "Robust Multipath-Assisted SLAM with Unknown Process Noise and Clutter Intensity" Remote Sensing 13, no. 9: 1625. https://doi.org/10.3390/rs13091625

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop