{'fraction_non_alphanumeric': 0.05603214865904751, 'fraction_numerical': 0.01191324422150654, 'mean_word_length': 4.107423580786026, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data Fusion for Multipath-Based SLAM: Combing Information from Multiple Propagation Paths

Multipath-based simultaneous localization and mapping (SLAM) is an emerging paradigm for accurate indoor localization with limited resources. The goal of multipath-based SLAM is to detect and localize radio reflective surfaces to support the estimation of time-varying positions of mobile agents. Radio reflective surfaces are typically represented by so-called virtual anchors (VAs), which are mirror images of base stations at the surfaces. In existing multipath-based SLAM methods, a VA is introduced for each propagation path, even if the goal is to map the reflective surfaces. The fact that not every reflective surface but every propagation path is modeled by a VA, complicates a consistent combination"fusion"of statistical information across multiple paths and base stations and thus limits the accuracy and mapping speed of existing multipath-based SLAM methods. In this paper, we introduce an improved statistical model and estimation method that enables data fusion for multipath-based SLAM by representing each surface by a single master virtual anchor (MVA). We further develop a particle-based sum-product algorithm (SPA) that performs probabilistic data association to compute marginal posterior distributions of MVA and agent positions efficiently. A key aspect of the proposed estimation method based on MVAs is to check the availability of single-bounce and double-bounce propagation paths at a specific agent position by means of ray-launching. The availability check is directly integrated into the statistical model by providing detection probabilities for probabilistic data association. Our numerical simulation results demonstrate significant improvements in estimation accuracy and mapping speed compared to state-of-the-art multipath-based SLAM methods.


I. INTRODUCTION
Emerging sensing technologies and innovative signal processing methods exploiting multipath propagation will lead to new capabilities for autonomous navigation, asset tracking, and situational awareness in future communication networks. Multipath-based simultaneous localization and mapping (SLAM) is a promising approach for obtaining position information of transmitters and receivers as well as information on their propagation environments. In multipath-based SLAM, specular reflections of radio signals at flat surfaces are modeled by virtual anchors (VAs) that are mirror images of base stations, called physical anchors (PAs) [1]. Typically, the number of VAs and their positions are unknown. Multipathbased SLAM methods can detect and localize VAs and jointly E. Leitinger  estimate the time-varying position of mobile agents [1]- [5]. The availability of VA location information makes it possible to leverage multiple propagation paths of radio signals for agent localization and can thus significantly improve localization accuracy and robustness [6]- [9].
A complicating factor in feature-based SLAM is measurement origin uncertainty, i.e., the unknown association of measurements with features [3]- [5], [23], [28]. In particular, (i) it is not known which map feature has generated by which measurement, (ii) there are missed detections due to low SNR or occlusion of features, and (iii) there are false positive measurements due to clutter. Thus, an important aspect of multipath-based SLAM is data association (DA) between these measurements and the VA. Probabilistic data association can increase the robustness and accuracy of multipath-based SLAM but introduces additional unknown parameters. Stateof-the-art methods for multipath-based SLAM are Bayesian estimators that perform the sum-product algorithm (SPA) on a factor graph [3]- [5] to avoid the curse of dimensionality related to the high-dimensional parameters space. Since the measurement models of multipath-based SLAM are nonlinear, most methods typically rely on sampling techniques [1]- [6]. Recently, multipath-based SLAM based on data collected by radios with ultra-wide bandwidth [29] or multiple antennas [2], [4], [30] was applied to indoor scenarios.
In existing methods for multipath-based SLAM, each VA represents a single propagation path from an agent to a PA. In particular, even if a reflective surface takes part in multiple propagation paths, each path is represented by a VA, and each VA is estimated independently [1]- [6], [31], [32]. The fact that not every reflective surface but every propagation path is modeled by a VA complicates a consistent combination "fusion" of statistical information across multiple propagation paths and PAs. It thus limits the accuracy and mapping speed of existing multipath-based SLAM.

B. Contributions and Notations
In this paper, we introduce a new statistical model that considers single-bounce and double-bounce propagation paths and enables data fusion for multipath-based SLAM. Contrary to existing multipath-based SLAM models, which represent propagation paths as SLAM features called VAs, in the proposed model, every reflective surface is represented by a SLAM feature referred to as master virtual anchor (MVA). Following the new MVA-based model, a factor graph is established, and an extension of the SPA multipath-based SLAM methods in [3]- [5] is developed. This SPA makes it possible to infer reflective surfaces by combining information provided by multiple propagation paths and PAs. We also present a particlebased implementation of the proposed SPA that performs raylaunching to determine the availability of each single-bounce and double-bounce propagation at a specific agent position. The availability check is directly integrated into the statistical model, i.e., it determines detection probabilities used for probabilistic data association. The resulting multipath-based SLAM method can provide fast and accurate estimates of a large number of reflective surfaces and mobile agent positions.
The key contributions of this paper are as follows.
• We introduce a statistical model and factor graph for multipath-based SLAM that considers single-and doublebounce propagation paths and facilitates data fusion.
• We extend the SPA for multipath-based SLAM [3]- [5] to the introduced factor graph to establish data fusion for multipath-based SLAM.
• We present a particle-based implementation that directly performs ray-launching to determine the availability of each single-bounce and double-bounce propagation path.
• We demonstrate significant improvements in estimation accuracy and mapping speed compared to existing multipath-based SLAM methods.
This paper advances over the preliminary account of our method provided in the conference publications [33], [34] by (i) extending the MVAs model to double-bounce propagation paths, (ii) introducing ray-launching to determine the availability of propagation paths, (iii) presenting a detailed derivation of the factor graph, (iii) discussing a particlebased implementation, and (iv) demonstrating performance advantages compared to the classical multipath-based SLAM [3], [4] in different scenarios.
Notation: Random variables are displayed in sans serif, upright fonts; their realizations in serif, italic fonts. Vectors and matrices are denoted by bold lowercase and uppercase letters, respectively. For example, a random variable and its realization are denoted by x and x, respectively, and a random vector and its realization by x and x, respectively. Furthermore, x and x T denote the Euclidean norm and the transpose of vector x, respectively; ∝ indicates equality up to a normalization factor; f (x) denotes the probability density function (PDF) of random vector x (this is a short notation for f x (x)); f (x|y) denotes the conditional PDF of random vector x conditioned on random vector y (this is a short notation for f x|z (x|z)). denoted as atan2(p 2 , p 1 ). The cardinality of a set X is denoted as |X |. Finally, δ(e) denotes the indicator function of the event e = 0 (i.e., δ(e) = 1 if e = 0 and 0 otherwise).

II. GEOMETRICAL RELATIONS
We consider a mobile agent equipped with an H-element antenna array and J PAs with known positions p where J is assumed to be known, in an environment described by S reflective surfaces with index s ∈ S {1, . . . , S}. At every discrete time step n, the element locations are denoted by p (h) n , h ∈ {1, . . . , H}, the agent position p n refers to the center of gravity of the array. We also define d (h) = p (h) n − p n and ψ (h) = atan2 p (h) 2,n − p 2,n , p (h) 1,n − p 1,n − ψ n , the distance from the reference location p n and the orientation, respectively, of the h-th element as shown in Fig. 1. At each discrete time slot n, the position p n ∈ R 2 and the array orientation ψ n of the mobile agent are unknown. The mobile agent transmits a radio signal, and the PAs act as receivers. 1 The radio signal arrives at the receiver via a line-of-sight (LOS) path as well as via multipath components (MPCs) originating from the reflection of surrounding objects.
We restrict the representation to MPCs related to singlebounce and double-bounce propagation paths. In particular, associated with PA j there are |D S | = S single-bounce path-related VAs [3], [35] at unknown positions p ss ,va ∈ R 2 with indexed (s, s ) ∈ D D {(s, s ) ∈ S × S|s = s }. Therefore, the maximum number of VAs for PA j is given by |D| = S + S(S − 1) with D {(s, s ) ∈ S × S} = D D ∪ D S . 2 By applying the image-source model [35], [36], a VA associated with a single-bounce propagation path is the mirror image of p (j) pa at reflective surface s ∈ S given as    where (s, s) ∈ D S , u s is the normal vector of reflective surface s, and e s is an arbitrary point on this surface. The second summand in (1) represents the normal vector w.r.t. to reflective surface s in direction u s with the length of two times the distance between PA j at position p pa . By applying the image-source model [35], [36] again on VA at position p where (s, s ) ∈ D D . For conveniently addressing PA-related variables and factors, we also define p  pa . An example is shown in Fig. 1. Note that the LOS path is associated with the PA itself. Due to blockage of propagation paths or geometrical restrictions between the reflective surfaces lead to availability constraints of propagation paths to certain VA positions [35]. Especially, double-bounce paths and their corresponding VA positions have limited availability depending on the mobile agent position p n (see Fig. 2b and Fig. 2c). Therefore, the number of available VAs for PA j is smaller than |D| = S+S(S−1). The proposed algorithm performs an availability check of each path using ray-launching [35] to determine whether a VA provides information as it is detailed in Section III-D.

A. MVA-based Model of the Environment
A reflective surface is involved in multiple propagation paths, thus giving rise to multiple VAs. To enable the consistent combination, i.e., "fusion" of map information provided by distance measurement of different PAs, we represent reflective surfaces by S unique MVAs at unknown positions p s,mva ∈ R 2 , s ∈ S. The unique MVA position p s,mva ∈ R 2 is defined as the mirror image of [0 0] T on the reflective surface S. By using some algebra, the transformation from MVA p s,mva to VA position p (j) ss,va can be obtained as where ·, · represents the inner-product between two vectors. Similarly, the inverse transformation from a VA to a MVA can be obtained Note that the inverse transformation in (4) is very useful to determine a proposal distribution for the MVA states as described in Section V. Details of the derivation are provided in the supplementary material. As an example, Fig. 2 shows three scenarios with two reflecting surfaces described by two MVAs at positions p 1,mva and p 2,mva . Fig. 2a shows a scenario with two PAs j ∈ {1, 2}, the according VAs, and an agent at position p n . Each PA generates one VA associated with a single-bounce propagation path. Fig. 2b shows a scenario with one PA, the according VAs associated with single-bounce and double-bounce propagation paths, and an agent at two positions p n1 and p n2 . Note that in scenarios with perpendicular surfaces, the order of the surfaces, i.e., surface s -surface s or surface s -surface s does not influence the calculation of the VA-position associated with the double-bounce propagation path. Depending on the agent position p n , only one double-bounce propagation path (related to one of the two orders) is available (Section III-D). Fig. 2b shows a scenario with non-perpendicular surfaces (acute angle between the surfaces). In this case, different order of the surfaces leads to difference VA-positions. Note that in such scenarios exist regions of the agent position p n (dependent on the PA positions and the angle between the surfaces) for which both double-bounce propagation path (associated with the VA-positions) are available (Section III-D).

III. SYSTEM MODEL
At each time n, the state of the agent x n = [p T n v T n ψ n ], where v n is the agent velocity. We assume that the array is rigidly coupled with the movement direction. As in [3], [37], [38], we account for the unknown number of MVAs by introducing PMVAs s ∈ S n {1, . . . , S n }. The number S n of PMVAs is the maximum possible number of actual MVA that produced a measurement so far [3], [38] (where S n increases with time). PMVA states are denoted as y s,n = p T s,mva r s,n T . The existence/nonexistence of PMVA s is modeled by the existence variable r s,n ∈ {0, 1} in the sense that PMVA s exists if and only if r s,n = 1. It is considered also if PMVA s is nonexistent, i.e., if r s,n = 0. The states p T s,mva of nonexistent PMVAs are obviously irrelevant. Therefore, all PDFs defined is an arbitrary "dummy PDF" and f s,n ∈ [0, 1] is a constant and can be interpreted as the probability of non-existence [3], [38].

A. Measurements and New PMVAs
The distance and AOA measurements 3 related to the path agent at position p n -VA at position p (j) ss ,va with (s, s ) ∈D n are given by where ν τ m,n p s,mva , p s ,mva , p n . Note that with the proposed MVA-based measurement model, at time n, the measurements collected by all PAs j ∈ {1, . . . , J} can provide information on the same MVA and agent positions p s,mva , s ∈ S and p n , respectively.
It is assumed that an MVA generates at most one measurement and that a measurement originates from at most one MVA. PA j at position p s,mva is related to the detection of a measurement within the pre-estimation stage as well as the availability check of a single-bounce propagation path between the agent position p n and the VA position p (j) ss,va with (s, s) ∈ D S,n . The same holds for the double-bounce propagation path between the agent position p n and the VA position p m,n may also not originate from any MVA. This type of measurement is referred to as a false positive. A Poisson point process models false positve measurements with mean µ fp and PDF f fp (z (j) m,n ). Newly detected MVAs, i.e., MVAs that generated a measurement for the first time, are modeled by a Poisson point process with mean µ n and PDF f n (p s,mva ). Newly detected MVAs are represented by new PMVA states y n } in our statistical model [3], [38]. Each new PMVA state corresponds to a measurement z (j) m,n ; r s,n = 1 implies that measurement z (j) m,n was generated by a newly detected MVA. This implies that each measurement z (j) m,n is assumed to originate from a single-bounce path corresponding to exactly one PMVA s ∈ D S,n (not a pair of PMVAs), even though some of the measurements are potentially related to double-bounce paths. However, new PMVA states y (j) n,k generated by measurements that are related to double-bounce path have a low probability of being associated again to measurements. Therefore their existence probability will drop immediately to very low values. We also introduce by y (j) n ,n T the joint vector of all new PMVA states related to measurements. Introducing new PMVA for each extracted measurement leads to several PMVA states that would grow with time n. Thus, to keep the proposed SLAM algorithm feasible, a sub-optimum pruning step is performed, removing PMVAs with a low probability of existence; this will be further discussed in Sec. IV-A.

B. Legacy PMVAs and State Transition
At time n, measurements are incorporated sequentially across PAs j ∈ {1, . . . , J}. Previously detected MVAs, i.e., MVAs that have been detected either at a previous time n < n or at the current time n but at a previous PA j < j, are represented by legacy PMVA states y (j) s,n . New PMVAs become legacy PMVAs when the next measurements-either of the next PA or at the next time instance-are taken into account. In particular, the MVA represented by the new MVA state y (j ) m ,n introduced due to measurement m of PA j at time n n is represented by the legacy PMVA state y (j) s,n at time +m . The number of legacy PMVA at time n, when the measurements of the next PA j are incorporated, is updated according to where S (1) n is equal to the number of all measurements collected up to time n and PA j−1. The vector of all legacy PMVA states at time n and up to PA j can now be written as y (j) Let us denote by y (1) n y T 1,n · · · y T Sn−1,n T , the vector of all legacy PMVA states before any measurements at time n have been incorporated. After the measurements of all PAs j ∈ {1, . . . , J} have been incorporated at time n, the total number of PMVA states is and the vector of all PMVA states at time n is given by y n = Note that this sequential incorporation of new PMVA states is based on the multisensor multitarget tracking approach introduced in [38, Sec. VIII]. We also define the number of VAs for each PA given as |D Legacy PMVAs states y s,n and the agent state x n are assumed to evolve independently across time according to state-transition PDFs f y s,n y s,n−1 and f (x n |x n−1 ), respectively. If PMVA k exists at time n −1, i.e., r s,n−1 = 1, it either disappears, i.e., r s,n = 0, or survives, i.e., r s,n = 1; in the latter case, it becomes a legacy PMVA at time n. The probability of survival is denoted by p s . Suppose the PMVA survives. In that case, its position remains unchanged, i.e., the state-transition pdf of the MVA positions p s,mva is given by f p s,mva p s,mva = δ p s,mva − p s,mva . Therefore, f p s,mva , r s,n p s,mva , r s,n−1 for r s,n−1 = 1 is obtained as If MVA s does not exist at time n−1, i.e., r s,n−1 = 0, it cannot exist as a legacy PMVA at time n either, thus we get f p s,mva , r s,n p s,mva , r s,n−1 = 0 = f D p s,mva , r s,n = 0 0, r s,n = 1.
For j 2, we also define f (j) y (j) s,n y (j−1) and It it assumed that at time n = 0 the initial prior PDF f y s,0 , s = 1, . . . , S 0 and f (x 0 ) are known. All (legacy and new) PMVA states up to time n are denoted as y 0:n y T 0 · · · y T n T and all agent states up to time n are denoted as x 0:n

C. Data Association Uncertainty
Mapping of reflective surfaces modeled by MVA is complicated by the data association uncertainty: at time n it is unknown which measurement z D,n associated with single-bounce and double-bounce propagation path. Any feature-to-measurement association (which considers associations to single PMVAs and MVA-MVA pairs as well as to PA j itself) is described by feature-oriented association variables generates measurement m 0, if legacy feature ss does not generate any measurement (13) with (s, s ) ∈D (j) n and stacked into the feature-oriented association vector as a n ,n T . To reduce computation complexity, following [3], [37]- [41], we use a redundant description of PMVA-measurement associations, i.e., we introduce measurement-oriented association variables if measurement m is not generated by any legacy feature ss (14) and stacked into the measurement-oriented association vector as a (j) Note that any data association event that can be expressed by both a joint feature-oriented association vector a (j) n and measurement-oriented association vector a (j) n is a valid event in the sense that an PMVA generates at most one measurement. A measurement is originated by at most one PMVA. This hybrid representation of data association makes it possible to develop scalable SPAs for simultaneous agents localization and MVA detection, and tracking [3], [37], [38], [41]. Finally, we also introduce the joint association vectors a n = a

D. Detection Probability and Availability of Propagation Paths
The availability check is performed by establishing a ray [35], [42] between the agent at position p n and a VA at position p pa is available 0, path is not available (16) with (0, 0), d,ss,n , path from agent at p n to VA at p (j) ss,va is available 0, path is not available (17) with (s, s) ∈ D (j) S,n as well as d,ss ,n , path from agent at p n to PMVA at p (j) ss ,va is available 0, path is not available (18) with (s, s ) ∈ D (j) D,n .

IV. PROBLEM FORMULATION AND PROPOSED METHOD
In this section, we formulate the detection and estimation problem of interest and present the joint posterior PDF and factor graph underlying the proposed SLAM method.

A. Detection of MVAs and State Estimation
At each time n, a channel estimator [21], [23]- [27] (preestimation stage) provides M We aim to estimate the agent state x n using all available measurements z 1:n from all PAs up to time n. In particular, we calculate an estimatex n by using the minimum meansquare error (MMSE) estimator [43,Ch. 4] x n x n f (x n |z 1:n )dx n .
The map of the environment is represented by reflective surfaces described by PMVAs. Therefore, the positions p s,mva of the detected PMVAs s ∈ {1, . . . , S n } must be estimated. This relies on the marginal posterior existence probabilities p(r s,n = 1|z 1:n ) = f (p s,mva , r s,n = 1|z 1:n )dp s,mva and the marginal posterior PDFs f (p s,mva |r s,n = 1, z 1:n ) = f (p s,mva , r s,n = 1|z 1:n )/p(r s,n = 1|z 1:n p k,mva p k,mva f (p k,mva |r s,n = 1, z 1:n ) dp k,mva . (20) The calculation of f (x n |z 1:n ), p(r s,n = 1|z), and f (p k,mva | r s,n = 1, z 1:n ) from the joint posterior f (y 0:n , x 0:n , a 1:n , a 1:n |z 1:n ) by direct marginalization is not feasible. By performing sequential message passing using the SPA rules [3], [37], [44], [45] on the factor graph in Fig. 3, approximations ("beliefs")f x n andf s y s,n of the marginal posterior pdfs f (x n |z 1:n ), p(r s,n = 1|z 1:n ), and f (p k,mva | r s,n = 1, z 1:n ) can be obtained in an efficient way for the agent state as well as all legacy and new PMVAs states s ∈ S n .

B. The Factor Graph
By using common assumptions [3], [38], and for fixed measurements z 1:n , it can be shown that the joint posterior PDF of y 0:n , x 0:n , a 1:n , and a 1:n , conditioned on z 1:n is given by (15) as shown on top of the page, where we introduced the functions q P x n , a for (s, s) ∈ D (j) S,n by In case the joint feature-oriented association vector a n and the measurement-oriented association vector a n do not describe the same association event, at least one indicator function in (15) is zero. Thus f (y 0:n , x 0:n , a 1:n , a 1:n |z 1:n ) is zero as well. A detailed derivation of this factor graph is provided in the supplementary material. The factor graph representing factorization (15) is shown in Fig. 3.
V. PROPOSED SUM-PRODUCT ALGORITHM Since our factor graph in Fig. 3 has cycles, we have to decide on a specific order of message computation [44], [46]. We choose this order according to the following rules: (i) messages are only sent forward in time 4 ; (ii) messages are only sent from PA j −1 to PA j, i.e., the measurements of PA are processed serial, thus, PA j−1 establishes new PMVAs that are acting as legacy PMVAs for PA j; (iii) iterative message passing is only performed for data association [3], [4], i.e., in particular, for the loops connecting different PMVAs, we only perform a single message passing iteration; and (iv) along an edge connecting an agent state variable node and a new PMVA state variable node, messages are only sent from the former to the latter. With these rules, the generic message passing equations of the SPA [44] yields the following operations at each time step. The corresponding messages are shown in Fig. 3.
We note that similarly to the "dummy pdfs" introduced in Section III, we consider BP messages ϕ y k,n also for the non-existing PF states, i.e., for r (j) k,n = 0. We define these messages by setting ϕ a (j) k,n , 0 = ϕ (j) k,n (note that these messages are not pdfs and thus are not required to integrate to 1).

A. Prediction Step
First, a prediction step is performed for the agent and all legacy PMVAs s ∈ S n−1 . Based on the SPA rule, the prediction message for the agent state is given by  f p s,mva , r s,n p s,mva , r s,n−1 ×f p s,mva , r s,n−1 dp s,mva , s ∈ {1, . . . , S n−1 }, where the beliefs of the mobile agent state, f (x n−1 ), and of the PMVA states,f p s,mva , r s,n−1 , were calculated at the preceding time n−1. Inserting (9) and (10) for f p s,mva , r s,n p s,mva , r s,n−1 = 1 and f p s,mva , r s,n p s,mva , r s,n−1 = 0 , respectively, we obtain α s p s,mva , 1 = p s δ p s,mva − p s,mva f s p s,mva , 1 dp s,mva , and α s p s,mva , 0 = α s,n f D p s,mva with α s,n = (1− p s f s p s,mva , 1 dp s,mva +f s,n−1 .
We note thatf s,n−1 f s p s,mva , 0 dp s,mva approximates the probability of non-existence of legacy PMVA s at the previous time step n − 1.

B. Sequential PA Update
At iteration j ∈ {1, . . . , J}, the following operations are calculated for all legacy and new PMVAs.
2) Checking the Availability of Propagation Paths: The proposed SPA algorithm performs an availability check for each propagation path using ray-launching [35], [42] to determine whether a VA can provide map information. First, the VA positions p pa is performed based on the fast line intersection algorithm [42]. In case the path between the agent position and a VA position or between two VA positions (double-bounce path) intersects with a reflective surface (e.g., it is blocked), the corresponding path is not available at the agent position p n . Hence, the corresponding VA cannot be associated with measurements at time n.

3) Measurement Evaluation for the LOS Path:
The message β 00 a (j) 00,n passed from the factor node q P x n , a s,n , a (j) ss,n , x n ; z (j) n × dx n dp (j) s,mva + 1 a (j) ss,n α (j) s,n .
For s = s , the message β ss a ss ,n , x n ; z (j) n × dx n dp (j) s,mva + 1 a Using the expression of q S p (j) m,mva , r 6) Iterative Data Association: Next, from β a m,n are calculated iteratively according to [38], [41] for ss ∈D 7) Measurement Update for the Agent: : The messages γ (j) 00 (x n ) send from the factor node q P x n , a The messages γ (j) ss (x n ) send from the factor node q S p (j) s,mva , 1, a (j) ss,n , x n ; z (j) n to the agent variable node are given for s = s by ss,n q S p (j) s,mva , 1, a (j) ss,n , x n ; z (j) n × α s p (j) s,mva , 1 dp (j) s,mva + η a (j) ss,n = 0 α (j) and from the factor node q D p (j) ss ,n , x n ; z (j) n to the agent variable node are given for s = s by ss ,n , x n ; z (j) n × dp (j) s,mva dp (j) s ,mva + η a ss ,n , x n ; z (j) n × dx n dp (j) s ,mva Based on these messages, the message sent to the next PA γ y (j) k,n is given by s,mva , 1 = q S p (j) m,mva , 1, a (j) m,n , x n ; z (j)

C. Belief Calculation
After the messages for all PAs, j = {1, . . . , J} are computed, the belieff (x n ) of the agent state can be calculated as normalized production of all incoming messages [44], i.e., A computationally feasible sequential particle-based message passing implementation can be obtained along the line with [3], [37], [45]. In particular, we adopted the approach in [45] using the "stacked state" comprising the agent state and the PMVA states. To avoid that the number of PMVA states grows indefinitely, PMVAs states with p(r (j) s,n = 1|z 1:n ) below a threshold p pr are removed from the state space ("pruned") after processing the measurements of each PA j.

VI. SIMULATION RESULTS
The performance of the proposed MVA-based SLAM algorithm is validated and compared with the multipath-based SLAM algorithm from [3] 5 using synthetic measurements.

A. Common simulation setup
The following setups and parameters are commonly used for synthetic and real measurements. The agent's state-transition pdf f (x n |x n−1 ), with x n = [p T n v T n ] T , is defined by a linear, near constant-velocity motion model [47,Sec. 6.3.2], i.e., x n = Ax n−1 + Bw n . Here, A ∈ R 4×4 and B ∈ R 4×2 are as defined in [47,Sec. 6.3.2] (with sampling period ∆T = 1s), and the driving process w n is iid across n, zero-mean, and Gaussian with covariance matrix σ 2 w I 2 , where I 2 denotes the 2 × 2 identity matrix and σ w denotes the acceleration noise standard deviation. For the sake of numerical stability, we introduced a small regularization noise to the PMVA state p k,mva at each time n, i.e., p k,mva = p k,mva + ω k , where ω k is iid across k, zero-mean, and Gaussian with covariance matrix σ 2 a I 2 and σ a = 10 −3 m. We performed 500 simulation runs using 50000 particles, each using the floor plan and agent trajectory shown in Fig. 4. In each simulation run, we generated noisy distance and AOA measurements stacked into the vector z (j) m,n according to (5) and (6) The performance of the different methods discussed is measured in terms of the root mean-square error (RMSE) of the agent position, as well as the optimal subpattern assignment (OSPA) error [48] of VAs and MVAs. The OSPA takes into account data association uncertainty and penalizes an incorrect number of estimated objects. We calculate the optimal subpattern assignment (OSPA) errors based on the Euclidean metric with cutoff parameter c = 5 m and order p = 1. The mean optimal subpattern assignment (MOSPA) errors, mean RMSEs of each unknown variable are obtained by averaging over all simulation runs.

B. Performance in Non-Rectangular Room
In this experiment, we compare the complete version of the proposed MVA-based SLAM algorithm (vis.), which includes the availability checks as introduced in Sec. III-D and Sec. V-B2, with a reduced variant, where we deactivate availability checks (no vis.), i.e., p d (p n , p  Fig. 5a and 5b of the algorithm variant with availability check converge faster and to a much smaller value than those of the variant without availability check. This is because in the scenario investigated, several VAs corresponding to the left as well as the lower walls are not available over large parts of the trajectory (this is the case especially for VAs of p (2) pa ). See also Fig. 2c, which provides a graphical explanation. Thus, the algorithm variant without availability check tends to deactivate the corresponding MVAs (i.e., strongly lower the probability of existence) as some of the corresponding VAs, which are expected to be detected with p d are not observed for significant amounts of time. The RMSE of the agent position are not strongly influenced by this deactivation since still sufficient positon-related information is provided by the two PAs and the remaining VAs.

C. Comparison of factorized SLAM and MVA-based SLAM
In this experiment, we compare the complete version of the proposed MVA-based SLAM algorithm (prop.) and multipathbased SLAM algorithm from [3]. We consider the indoor scenario shown in Figure 6. For ease of comparison, we chose the scenario to be identical to [3]. The scenario consists of four reflective surfaces, i.e., K = 4 MVAs, as well as two PAs at positions p  For the proposed algorithm, only one of the 500 simulation runs diverged, but 4 % of the simulation runs diverged for the algorithm in [3]. As an example, Fig. 6 depicts for one simulation run the posterior PDFs represented by particles of the MVA positions and corresponding reflective surfaces as well as estimated agent tracks. We also present results where we vary the measurement noise standard deviations. To this end, we introduce factors f std, all and f std,los . Both factors equally increase all measurement standard deviations σ τ mean agent RMSE, respectively, over all time steps n as a function of f std,all or f std,los . In particular, in Fig. 8a, 8c and 8e we varied f std,all and kept f std,los 1 fixed, while in Fig. 8b, 8d and 8f we varied f std,los and kept f std,all 4 fixed. The MOSPA error of the proposed algorithm related to VAs and PAs is shown in Fig. 7a. It can be seen that the MOSPA drops significantly after only a few time steps. In contrast, the algorithm presented in [3] only converges rather slowly to a small mapping error. Furthermore, the MOSPA error of the proposed algorithm converges along the agent track to a smaller value, i.e., to a smaller mapping error. This shows that the proposed MVA-based SLAM algorithm efficiently exploits all measurements for the map features provided by the PAs. Fig. 7b shows the MOSPA error of the MVA positions, which confirms the results seen in Fig. 7a. Finally, the RMSE of the agent position in 7c is slightly increased. At the same time, the map is unstable (i.e., the MOSPA errors are high) but converges to a small value, only increasing slightly during changes of the agent's direction. In contrast, the RMSE of the algorithm from [3] even increases during the straight walk of the agent between n = 50 and n = 100. Both methods show similar performance in terms of the RMSE of the agent orientation. The VAs MOSPA errors in Figs. 8a to 8d emphasize these observations. Naturally, all error values increase with increasing f std,all . However, the MOSPA errors of the proposed method increase much slower due to the additional information contained in the second-order VAs, which are leveraged by the MVA model. In contrast, the agent RMSE in Fig. 8e remains constant for both methods as there is still enough information available for proper localization, which is mainly provided by the PAs. Thus, when additionally increasing f std,los in Fig. 8f, the RMSE of the algorithm from [3] significantly increases, while the agent RMSE of the proposed algorithm remains approximately constant due to the increased map stability.

VII. CONCLUSIONS
In this paper, we introduced data fusion for multipathbased SLAM. A key novelty of our approach is to represent each reflective surface in the propagation environment by a single MVA. This representation addresses a key limitation of existing multipath-based SLAM methods that represent every propagation path by a VA which complicates a consistent combination of statistical information across multiple paths and base stations and, in turn, limits the accuracy and mapping speed of existing multipath-based SLAM methods. A key aspect in leveraging the advantages of the introduced MVAbased model, was to check the availability of single-bounce and double-bounce propagation paths at a specific agent position by means of ray-launching. This availability check was directly integrated into the statistical model by providing detection probabilities for probabilistic data association. Our