=Paper=
{{Paper
|id=Vol-3248/paper31
|storemode=property
|title=Unsimultaneous Time of Arrival Multipath-Based Localization and Mapping
|pdfUrl=https://ceur-ws.org/Vol-3248/paper31.pdf
|volume=Vol-3248
|authors=Conghua Chen,Ao Peng,Xuemin Hong
|dblpUrl=https://dblp.org/rec/conf/ipin/ChenPH22
}}
==Unsimultaneous Time of Arrival Multipath-Based Localization and Mapping==
https://ceur-ws.org/Vol-3248/paper31.pdf
Unsimultaneous Time of Arrival Multipath-Based
Localization and Mapping
Conghua Chen1 , Ao Peng*1 and Xuemin Hong1
1
School of Informatics, Xiamen University, Xiamen, China
Abstract
Indoor positioning based on time of arrival (TOA) can be a huge challenge. The complexity and uncer-
tainty of data association (DA) due to the multipath effect and the visibility of anchor observation due
to continuous terminal movement are important issues to be addressed by current indoor positioning
algorithms. Multipath-assisted localization treats the reflected signal received at the base station as
the direct signal received at the virtual anchor point (VA), which can significantly improve the mul-
tipath interference problem. Unlike traditional multipath SLAM, anchors in the scenario are firstly
estimated in a training set when the terminal locations are known using Feature mapping single cluster-
probabilistic hypothesis density (FMSC-PHD) filtering. Then the terminal localization are solved with
a factor graph based belief propagation (BP) algorithm based on the estimated anchors. Experimental
results demonstrate the excellent performance of the algorithm in mapping and localization.
Keywords
Indoor positionings, TOA, PHD, RFS
1. Introduction
With the repaid development of mobile equipment and 5G networks, indoor positioning and
location-based services provide us with more convenience and raise many new requirements.
The applications for indoor positioning are very promising, including robotics[1], Internet of
Things[2], location-aware communication[3] and so on, all have a large demand. The indoor
environment has severe shading, multipath effect, and Doppler effect, which severely limit the
accuracy and reliability of positioning. In recent years, many solutions have also been derived
including inertial guidance, geomagnetic, LiDAR, radio, and other measurement sources that
can be used for localization[4]. Among them, radio signals have become the focus of research
with their advantages of easy deployment, wide applicability, and low cost.
Wireless signals like Bluetooth, UWB, Wi-Fi, and other positioning information source are
widely available in our indoor space[5]. Currently, there are two major radio-based localization
approaches[2, 6, 7]: fingerprint-based and geometry-based methods. Fingerprinting localization
is a method of correlating the specific location with signal features[8]. On the other hand, the
geometry-based approach has great potential for generalization. Geometric methods consist
of range-based and angle-based methods. TOA based on the time of arrival and AoA based
on the angle of arrival is commonly used measurements respectively. Although it can solve
IPIN 2022 WiP Proceedings, September 5 - 7, 2022, Beijing, China
$ chenconghua11@gmail.com (C. Chen); pa@xmu.edu.cn (A. Peng*); xuemin.hong@xmu.edu.c (X. Hong)
© 2022 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR
Workshop
Proceedings
http://ceur-ws.org
ISSN 1613-0073
CEUR Workshop Proceedings (CEUR-WS.org)
the problems of spatial consistency in indoor positioning, it requires high angular resolution
of terminals[8]. This means that to achieve better accuracy, the complexity and cost of the
equipment will be greatly increased. In contrast, TOA is a very accessible measurement that
requires only very simple equipment to obtain observations of considerable accuracy. TOA is
currently the most likely option for large-scale deployment.
Multipath-assisted positioning[9, 10] is one of the best solutions for indoor positioning to
solve the multipath effect. It assumes that the walls are all smooth, so the multipath components
are the specular reflection. Specular reflection components are generated in the process of
transmitting radio signals from the mobile terminal to the base station, and they are equivalent to
the Line of Sight(LOS) signals generated by the VAs. Multipath-assisted localization transforms
the previous multipath interference into LOS path signals sent by multiple virtual base stations.
Multipath-assisted localization takes full advantage of the geometric properties of space, not only
overcoming multipath interference but also adding many useful signal sources for localization.
Although multipath-assist localization can take full advantage of the otherwise interfering
multipath reflection components, there are still some problems with the accuracy and complexity
of DA and visibility of anchors in the region of interest. Since there is no distinguishability
among multipath components, we cannot know the relation between anchors and multipath
components. TOA-based indoor localization requires anchors’ position and their corresponding
measurements to solve the location of the terminals. Therefore, a method to solve DA is needed.
In addition, the indoor environments are often heavily obscured, which can cause visibility
problems the terminal may not observe the PAs and VAs when it moves.
Our approach is divided into two phases, first estimating the anchor in the environment using
FMSC-PHD, and then localizing the terminal when the anchor prior information is available.
Our method split into two phases since it is difficult for SC-PHD to perform SLAM with only
TOA observations[11]. When SLAM is performed by SC-PHD, the estimated localization of the
terminal depends on being able to estimate the approximate location of the anchors in a short
time. Estimated the accurate position of anchors needs a process to convergence when SC-PHD
is performed on SLAM with only TOA metric. This means that neither the terminals nor the
anchors can be reliably located. Even though [12, 13] solved this problem to some extent, it is
still limited by the observation of the LOS component. Therefore, by adding terminal locations
as a priori information in the training set, FMSC-PHD can estimate the location of anchors in
the scenario using the terminal locations and TOA.
FMSC-PHD is a filter method extended from PHD[14] and SC-PHD[15, 16] which directly
avoids the complex DA process and solves the uncertainty problem of it. In addition, FMSC-PHD
uses the RFS for modeling anchors in the scenario, which solves the problem of the visibility
of the anchor during the motion of the terminal. After estimating the anchors’ position in the
training set, the localization of the terminal in this indoor scenario can be reduced to the DA
between TOA observations and estimated anchors and the estimation of the terminal positioning
based on the mapping result. We build a probabilistic model of DA and terminal positioning
based on factor graphs and solve it with the BP algorithm. The BP algorithm running on factor
graphs has a great advantage for solving edge probabilities and is suitable for solving DA and
location estimation.
Key innovative contributions of this paper include the following:
VA
PA
Terminal
Figure 1: Example of PA and VA
• We divide the complex SLAM process into two independent processes, first using FMSC-
PHD to obtain the features of the scenario such as PAs/VAs, and then performing terminal
positioning based on these estimated anchors. The limitations of SC-PHD in TOA mea-
surements are avoided, and the complexity of the system is reduced at the same time.
• We build a Bayesian model based on FMSC-PHD filtering to represent PA/VA in the form
of RFS, which well solves the problem of visibility of anchor during terminal moving. the
FMSC-PHD inherited from PHD filtering also avoids the process of DA
• Our proposed algorithm is verified through simulation experiments, and the experimental
results demonstrate the high accuracy and scalability of our algorithm.
2. Problem Formulation
As shown in Fig.1, we consider a static PA as a receiver and a mobile terminal as a transmitter
in a two-dimensional space (2-D). The terminal is constantly transmitting signals to the base
station during the movement. The base station is also receiving multipath signals to extract TOA
parameters. In multipath-assisted localization, a PA in a scenario may correspond to multiple
mirrored VAs. From the principle of geometric optics, it is known that in the premise of specular
reflection, the reverse extension of the reflected signal intersects with the wall-normal through
PA at the same point, which is the mirror point of PA, and the linear distance between the
terminal and the mirror point is equal to the reflected signal propagation distance. We only
consider first-order reflections, and subsequently no longer make a strict distinction between
PA and VAs.
The state of the mobile terminal at step k is represented by a vector xk = [pk , vk ], where
pk = [xk,1 , xk,2 ] denotes position and vk = [vk,1 , vk,2 ] denotes speed. The position of the nth
static anchor at step k can be expressed as mk,n = [mk,1 , mk,n,2 ]. During the movement of the
terminal, encountering obstacles and long distances may lead to missed detections, and new
anchors may be detected at the same time, which causes visibility problems. The number of
anchors which can be observed by the terminal is variable. RFS is suitable for representing such
a variable number of observed anchors. Hence all anchor nodes observed at each moment can
be represented as an RFS, with each set element representing a random variable of the anchor
distribution. The cardinality of RFS represents the number of anchors at the moment. The RFS
containing all environmental features at step k can then be expressed as
{︀ }︀
Mk = Mk−1 ∪ Bk = mk,1 , . . . , mk,|Mk | , (1)
where Mk−1 = {mk,1 , ..., mk,nk } is anchor set survive at step k-1, Bk represents newly
detected anchors set. Bk can be modeled as a Possion point process with intensity bk = µk I (·)
where I (·) is the distribution of new birth anchors and the new birth rate µk means the average
number of new birth anchors.
Assuming terminal moving with constant velocity obeys a linear Gaussian model. Its state
equation can be express as
[︂ ]︂ [︂ ∆T 2 ]︂
I2 ∆T · I2 2 · I2
xk = Fxk−1 + Guk−1 = xk−1 + uk−1 , (2)
02 I2 ∆T · I2
where uk−1 is accelerate noise, IN and 0N represents N dimension identity matrix and zero
matrix respectively, ∆T represents the sample period.
Assume that an anchor corresponds to only one multipath component. Since wireless signals
are multipath propagated in indoor space, we can solve the parameters corresponding to different
multipath signals from the received signals. Provided that clocks between the terminal and base
station are synchronized, the signal s (t) received by the base station from the terminal r (t)
can be expressed as
|Mk |
∑︁
rk (t) = wn,k s (t − τn,k ) + dk (t) + nAW GN (t). (3)
n=1
The first term on the right-hand side of the equation represents |Mk | multipath component
at step k. It can be considered the signal received by different anchors at step k. wn,k and
τn,k represents the nth multipath component complex amplitude and latency respectively. In
the ideal case the time delay is in proportional to the distance, which can be expressed as
τn,k = ‖xk − mk,n ‖ /c, where c is the speed of light. The second term dk (t) represents the
scattering components which affect the observation of specular reflections in the form of false
alarms. The last term nAW GN (t) is additive Gaussian white noise(AWGN). Based on the signal
model, the measurement nth zk,n between terminal and anchors can define as
√︁
zk,n = (xk,1 − mk,n,1 )2 + (xk,2 − mk,n,2 )2 + nk,n , (4)
where nk,n is the observation noise. Since the number of observations at step k is in dependent
on the detected anchors and false alarms. Hence the RFS of observation can be defined as
Zk = Z(xk , Mk ) ∪ Ck . (5)
The observations at step k consist of two sets. One is the TOA-based observation of the distance
with added noise. The miss-detection is mainly represented by the number of sets that may be
less than the number of current anchor points. The second represents the set Ck of false alarm
measurements caused by the scattering component. Ck can be considered as a Poisson point
process with intensity c(z) = λc U (z), where λc is the false alarm rate (also understood as the
average number of false alarms), and U (z) denotes a uniform distribution of false alarms over
the detection range.
The anchors estimated in the training set can assist the mobile terminal positioning. State
equation and measurement equation identical to those of the previous training set. The problem
of terminal localization then reduce to DA and location estimation. DA is described by two
vectors - a feature-oriented vector ak,i and an observation-oriented vector bk,m .
{︃
m ∈ ℳk , ith anchor generate mth measurement
ak,i =
0 , miss detect ith anchor
{︃
i ∈ ℐk , mth anchor generate ith measurement
bk,m =
0 , mth measurement is false alarm
where ℳk and ℐk correspond to measurement set and anchor set, respectively.
After obtaining the probability distribution of the DA, the prediction and update formulas of
the Bayesian filter can be used to recursively work out the distribution of the real-time terminal
location distribution. The Bayesian filter is updated in real-time, so the terminal positions can
be estimated in real time.
3. Mapping scenario Via FMSC-PHD
We first estimate the anchors in the scenario using FMSC-PHD. In the training set, the state
of the terminal is known, which can simplify the otherwise more complex SC-PHD that is
more applicable to our problem formulation. We also give the Monte Carlo implementation for
FMSC-PHD.
3.1. FMSC-PHD
According to SC-PHD[15], the prediction of terminal and anchors can be written as
∫︁
(︀ )︀ (︀ )︀
Dk|k−1 xk|k−1 , mk|k−1 = D̃k|k−1 mk|k−1 |xk−1
(︀ )︀ (6)
× sk−1 (xk−1 ) φk|k−1 xk|k−1 |xk−1 dxk−1
(︀ )︀
where sk−1 (xk−1 ) is terminal probabilistic distribution at step k-1 and φk|k−1 xk|k−1 |xk−1
(︀ )︀
is terminal state dynamics follow Markov process. D̃k|k−1 mk|k−1 |xk−1 is anchors’ PHD
conditioned on terminal state at step k.
The terminal state of the training set is known so that FMSC-PHD filter prediction can be
generated from (6). With the help of the sample property of the Dirac delta function, terminal
probability distribution s (xk ) also can be used for the case where the terminal state is known,
i.e. (︁ )︁
′
sk (xk ) = δ xk − xk (7)
where δ (·) is Dirac delta function. Substituting (7) into (6) and according to sampling property
of Dirac delta function, the joint PHD of terminal and anchors can be write as
(︀ )︀ (︁ ′
)︁ (︁ ′
)︁
Dk|k−1 xk|k−1 , mk|k−1 = φk|k−1 xk|k−1 |xk−1 D̃k|k−1 mk|k−1 |xk−1 (8)
According to (8), it can be seen that the joint FMSC-PHD of terminal and anchors predictions
at step k is the state transition equation multiplied by the prediction of conditional PHD of
anchors. When the state of the terminal is known, the terminal state at both step k-1 and
step k is a fixed value, and both can be described by a Dirac delta function like (7). Because
the uncertainty in the terminal locations is eliminated, the state transition equation does not
impact the anchor’s PHD prediction. Therefore, this joint distribution only needs to consider
the conditional PHD of anchors in a scenario where the terminal state is known. Hence the
FMSC-PHD prediction can be represented as
(︀ )︀ (︁ ′
)︁ (︁ ′
)︁
Dk|k−1 xk|k−1 , mk|k−1 = δ xk|k−1 − xk|k−1 D̃k|k−1 mk|k−1 |xk−1 . (9)
With previous processing, the prediction of FMSC-PHD was transformed into the prediction
of PHD(︀for anchor conditioned
)︀ on the terminal state alone. The prediction PHD of anchors
D̃k|k−1 mk|k−1 |xk−1 in step k can be specifically be expressed as[15]
(︀ )︀ (︀ )︀
D̃k|k−1 mk|k−1 |xk−1 = γk|k−1 mk|k−1 |xk−1
∫︁
(︀ )︀
+ D̃k−1 (mk−1 |xk−1 ) pS (mk−1 |xk−1 ) × φk|k−1 mk|k−1 |mk−1 ; xk−1 dmk−1
(10)
(︀ )︀
It consists of two parts, the new birth anchors γk|k−1 mk|k−1 |xk−1 and the prediction of the
(︀ )︀
anchor survive from the previous step. D̃k−1 (mk−1 |xk−1 ) and φk|k−1 mk|k−1 |mk−1 ; xk−1
represents the updated PHD of anchor at step k-1 and Markov transition probability of anchor.
The joint update function of SC-PHD can be shown as
sk|k−1 (xk ) LZk (xk )
Dk|k (xk , mk ) = ∫︀ D̃ (mk |xk ) (11)
sk|k−1 (xk ) LZk (xk ) dx k|k
where sk|k−1 (xk ) is the predicted terminal distribution, LZk (xk ) is the measurement likelihood
function and D̃k|k (mk |xk ) is the updated PHD of anchor conditional on terminals
The terminal part of the update formula can be simplified using known terminal positions.
The predicted terminal state is the same as the terminal position at step k. Hence, the distribution
of predicted terminal states can be replaced by a known Dirac delta function of the terminal
state. Substituted predicted terminal state Dirac delta function into (11), we can get
δ (x − x′ ) L (x′ )
(︀ k ′ )︀k Zk D̃k|k mk |x′k = δ xk − x′k D̃k|k mk |x′k
(︀ )︀ (︀ )︀ (︀ )︀
Dk|k (xk , mk ) = ∫︀
δ xk − xk LZk (xk ) dxk
(12)
Dirac delta function changes the integral of the denominator into a product, while the
likelihood function becomes the value of the function at the terminal position at step k. Similarly,
the likelihood function of the numerator is affected by the Dirac delta function and becomes
the value of the terminal state as a variable at step k. At this point, the values of the likelihood
functions of the numerator and denominator can be eliminated from each other, leaving the
updated part about the anchor PHD[15]
⎡ ⎤
∑︁ g (z|mk , xk ) pD (mk |xk )
D̃k|k (mk |xk ) = D̃k|k−1 (mk |xk ) × ⎣(1 − pD (mk |xk )) + ⎦,
ηz (mk |xk )
z∈Zk
∫︁ (13)
ηz (xk ) = κk (z) + D̃k|k−1 (mk |xk ) pD (mk |xk ) g (z|mk , xk ) dmk (14)
where pD (mk |xk ) represents the detection probability, g (z|mk , xk ) represents the likelihood
function, and κk (z) represents the false alarm. Likelihood factor in (13) has two parts: unde-
tected anchors and detected anchors based on measurements RFS Zk .
The advantage of PHD is that it does not require DA. We can see DA hidden in the likelihood
function, which can work out anchors’ PHD without calculating DA. With the training set
terminal position known, the PHD of the anchors at step k can be computed recursively at step
k-1 by a prediction and an update. The prediction equation and update equation reveal that the
FMSC-PHD constructed by the cluster model is very similar to the PHD filter[14] under the
condition that the terminal state is known. The difference is that the previous PHD was used to
estimate targets, but now we can use it to estimate anchors in the surrounding.
3.2. Monte Carlo implements of FMSC-PHD
Monte Carlo particle filter can be use to implements FMSC-PHD we performed in 3.1. The
updated PHD at step k-1 of anchors Dk−1 (mk−1 ) can be simulated by set of weighted particles
(i) (i)
{wk−1 , mk−1 }L k
i=1
Lk−1
∑︁ (i) (︁ (i)
)︁
Dk−1 (mk−1 ) ≈ wk−1 δ mk−1 − mk−1 , (15)
i=1
(i)
where Lk−1 is the number of particles and wk−1 is the weight of the ith particle at step k-1
respectively. Substituting particles simulation (15) into (10), we can get particles of predicted
PHD of anchors at step k
Lk−1 (︁ )︁ (︁ )︁
(i) (i) (i)
∑︁
Dk|k−1 (mk−1|k ) ≈ Dγ (mk ) + wk−1 ps,k mk−1 φk|k−1 mk|k−1 |mk−1 (16)
i=1
(︁ )︁ (︁ )︁
(i) (i)
where ps,k mk−1 is the probability that particles at moment k-1 survives at step k, φk|k−1 mk|k−1 |mk−1
is the transition probability of the ith particle from step k-1 to the step k. The predicted PHD of
anchors can be calculated by importance sampling.
Lk−1 (︁ )︁ φ (︁ )︁ D (m )
(i) (i) k|k−1 (i) γ k
(︀ )︀ ∑︁
Dk|k−1 mk|k−1 ≈ wk−1 ps,k mk−1 ×qk mk |mk−1 , Zk + qγ,k (mk |Zk )
qk qγ,k
i=1
(17)
Hence, we can get the approximate predicted PHD of anchors.
Lk−1 +Lγ (︁ )︁
(i) (i)
(︀ )︀ ∑︁
Dk|k−1 mk|k−1 ≈ wk|k−1 δ mk|k−1 − mk|k−1 (18)
i=1
where ⎧ (︁ )︁
⎨q m |m(i) , Z , i = 1, · · · , Lk−1 ,
(i) k k k−1 k
mk|k−1 ∼ (19)
⎩q (m |Z ) , i = Lk−1 + 1 · · · Lk−1 + Lγ,k ,
k k k
⎧
(i)
⎪
⎪ φk|k−1 ps,k wk−1
⎪
⎪ (︁ )︁ , i = 1, · · · , Lk−1 ,
⎨ qk mk |m(i) , Zk
⎪
⎪
(i) k−1
wk|k−1 ∼ (︁ )︁ (20)
(i)
D m
⎪
⎪
⎪
⎪ γ,k k
⎪
⎪ , i = Lk−1 + 1 · · · Lk−1 + Lγ,k .
Lγ,k qk (mk |Zk )
⎩
Substituted the predicted PHD particles set into (13) can get particle simulated updated PHD
Lk−1 +Lγ (︁ )︁
(i) (i)
(︀ )︀ ∑︁
Dk|k mk|k ≈ wk|k δ mk|k − mk|k . (21)
i=1
The anchor updated PHD particles are directly inherited from the surviving and newborn anchor
particles in the predicted PHD.
(i) (i)
mk|k = mk|k−1 . (22)
The ith particle weight is
(︁ )︁ (︁ )︁
(i) (i)
(︁ (︁ g z|mk|k−1 , xk pD mk|k−1 |xk
)︁)︁
(i) (i) (i) (i)
∑︁
wk|k = wk|k−1 1 − pD mk|k−1 + wk|k−1 (︁ )︁ .
(i)
z∈Zk η z m |x
k|k−1 k
(23)
It is possible to calculate the anchor PHD using the Monte Carlo method.
Since the anchors’ PHD are simulated by Monte Carlo methods using particles. Clustering
algorithm[17] can be taken to estimate anchors’ position from particles. The introduction of
the clustering algorithm brings instability and error to the resulting map, so it is crucial to find
a good estimation from all the training sets. However, we do not know the exact location of
the anchor nodes beforehand, so it is impossible to find the best-performing point directly by a
method.
To solve this tricky problem, we choose to use OSPA[18] as a measure for analyzing the
error of the anchor points. To overcome the unknown real anchors’ position, we make a small
change in its input. Since terminal trajectories are known, the distance can be calculated using
the terminal’s position and the anchor node’s estimated position at each moment. We can use
this distance and the observed distance as the OSPA input, and the result thus calculated can
be used as an alternative to the OSPA of estimated result and the exact anchors. We name it
measurement-oriented OSPA (MOSPA) for convenience.
4. localization based on Mapping result
With the help of the training set, we estimate the approximate locations of the anchor in the
scenario by FMSC-PHD filtering. The estimated anchors in the scenario can simplify the indoor
localization problem, which is then divided into two significant problems: DA and location
estimation of terminals. Solving the DA to obtain the correspondence between the observation
and the anchors enables further estimation of the anchor’s location by MMSE.
The relationship between observations and anchors in terms of probabilities usually involves
many operations for solving marginal probabilities. The direct calculation of edge probabilities
leads to the excessive complexity of the algorithm. We use factor graphs to represent the
probabilistic model of DA and then optimize the marginal probability solution in DA using
the BP algorithm. Since the localization of terminals is also based on the probabilistic form of
Bayesian filtering, it can also be represented by factor graphs. DA and terminal state distribution
can be computed directly and efficiently using the BP algorithm.
As mentioned earlier, DA between anchors and observations is described by feature-oriented
vector ak,i and observation-oriented vector bk,m . The prior state for Bayesian filtering at step k
can be expressed as p (xk , ak ). Through straightforward and well-known manipulation, the
likelihood function of DA and measurements conditional on the terminal state can be expressed
as[19]
|Mk |
(︃ (︀ ak,i )︀ )︃θd (ak,i ) |Zk |
p D (m |x
k,i k ) p z
(1 − pD (mk,i ))1−θd (ak,i ) ×
∏︁ ∏︁
p (Zk , ak |xk ) ∝ (︀ )︀ k ψc (ak,i , bk,j ),
i=1
κk zk,ak,i j=1
(24)
where |Mk | is the number of anchors in the region of interest, |Zk | is (︀the number )︀ of measure-
ments, pD (mk,i ) is the detect probability of ith anchor at step k, κk zk,ak,i represent false
alarm. θd (ak,i ) indicate whether the ith multipath component is detected or not, which can be
represent as {︃
i 0, aik = 0
θd (ak ) = (25)
1, aik ̸= 0.
Using Bayes’ rule and independence assumptions related to the prior probability density
function(pdf) and likelihood function, the joint posterior pdf of xk and ak at step k is obtained
Figure 2: Factor graph of (26)
as
|Mk | |Zk |
∏︁ ∏︁
p (xk , ak |Zk , mk ) ∝ ψp (xk , ak,i ) ψc (ak,i , bk,j ) (26)
i=1 j=1
where ⎧
⎪
⎨ [1 − pD (mk,i )] p (xk |Zk ) , ak,i = 0
(︀ )︀
ψp (xk , ak,i ) = pD (mk,i ) p zk,ak,i |xk , mk,i p (xk |Zk ) (27)
⎪ (︀ )︀ , ak,i ̸= 0
κk zk,ak,i
⎩
Since DA is based on the mutual constraints of the feature-oriented vector and the observation-
oriented vector, it can also be reformulated as
{︃
0, ak,i = j, bk,j ̸= i or bk,j = i, ak,i ̸= j
ψci,j (ak,i , bk,j ) = (28)
1, otherwise
From equation (26) we can obtain a factor graph as shown in Fig.2, and in turn we can run
the message propagation algorithm on the factor graph. The message propagated between ak,i
and bk,j can be obtained as
∑︁ ∏︁
µai →bj (bj ) = ψi (ai ) ψcij (ai , bj ) µbj ′ →ai (ai ) (29)
ai j ′ ̸=j
∑︁ ∏︁
µbj →ai (ai ) = ψcij (ai , bj ) µai′ →bj (bj ) (30)
bj i′ ̸=i
Since the factor graph has loops, there are no closed-form solutions. Approximate marginal
pdf can be obtained using iterative operations, and convergence was proved in the article[20].
The marginal probability density of the variable nodes can be computed by the computed
messages. The marginal pdf is associated with the location of the terminal, and the data at the
moment k can be expressed as
|Mk | ∫︁ |Zk |
∏︁ ∏︁
p (xk ) = ψp (xk , ak,i ) µbj →ai (ak,i )dak,i (31)
i=1 j=1
∫︁ |Zk |
∏︁
p (ak,i ) = ψp (xk , ak,i ) µbj →ai (ai ) dxk (32)
j=1
For estimating xk , we will develop an approximate calculation of the minimum mean-square
error (MMSE) estimator ∫︁
x̂MMSE
k = xk p (xk ) dxk (33)
5. experimental and Simulation result
In this section, to analyze the performance of the proposed FMSC-PHD filter and BP localization
algorithm, we apply it to simulation data within 2-D scenarios in Fig.1. The first situation
confiders the training set using the FMSC-PHD filter to position anchors. The second situation
considers the test set using the BP algorithm to simultaneously solve the DA and terminal
localization.
5.1. Analysis setup
State-Evolution Model The terminal’s state-transition pdf shown in Section 2 is defined
by a linear, near constant-velocity motion model[21] with sampling period ∆T = 1s. The
driving process uk is iid across k, zero-mean, and Gaussian with σu2 I2 accelerate noise, σu is
the accelerate noise. The anchors are static. However, implementing the FMSC-PHD algorithm
introduced a tiny driving process in the anchor state-evolution model for measurement noise.
Accordingly the state evolution is modeled as mk,n = mk−1,n + ωk,n , where ωk,n is iid across
k and n, zero-mean, and Gaussian with covariance matrix σm 2 I
2
Measurement Model According to the signal and measurement distance in Section 2, the
measurement noise nk,n is iid across k and n, zero-mean, and Gaussian with variance σz2 I2 . The
measurement model determines the likelihood function factors in (13) and (24).
Common Simulation Parameters The simulating teaching building, terminal trajectory
of the training set and test set, and static anchors show in Fig.1. The following parameters
are used for both the training set and test set. The false alarm measurements in (5) range
uniformly distribute on the region of interest, and the number follows the Poisson distribution.
The measurement we detect can be regarded as newly born anchors set in (1). The detection
probability is a constant value for anchors, i.e. pD (mk |xk ) = pD .
10
OSPA
OSPA map error (m)
8 MOSPA
6
4
2
0
0 50 100 150 200 250 300 350 400
10
8
no. of detected PF
6
4
Estimated Anchor
2
No. of Measurement
0
0 50 100 150 200 250 300 350 400
trajectory steps
Figure 3: The top figure plot training set results of OSPA and MOSPA of map error. The bottom figure
plot the number of estimated anchor and the number of measurement for training set.
Table 1
The 1 − σ error of cdf
Tests Benchmark Known Anchor Test1 Test2 Test3
Error 0.0134 0.0902 0.1469 0.4434
5.2. Result for Training Set
We used the common simulation parameters described above for our simulations based on a
training set. Terminal travels the diamond trajectory in Fig.4 to map the five anchors in the
scenarios with FMSC-PHD. We used pd = 0.9 and measurement standard variance σz = 0.1
and λc = 1 as parameters to analyze the FMSC-PHD algorithm. The PHD of newborn anchors
were each represented by 10000 particles.
Fig.3 shows the OSPA and MOSPA map error. The OSPA errors are based on the Euclidean
metric and use the cutoff parameter 10m and order 1. Fig.3 shows the simulated MOSPA can
approximate the real OSPA very well. These results demonstrate that the FMSC-PHD algorithm
can extract the best OSPA performance estimation in the scenarios. The OSPA meets a cutoff in
the 100 steps, corresponding to a big turn that solves the problem of space consist. We go on a
diamond trajectory, not a simple line. The estimated anchors are shown in Fig.4, which are very
close to the real anchor positions.
Fig.3 shows the number of anchors estimated and measurements generated. Due to visibility
problems, the clutter measurements and false alarms make the number of measurements vibrate
and offer larger than the actual number of anchors. Nevertheless, the performance of our
algorithm is still good. The number of estimated anchors is much more stable at around five,
the number of anchors shown in Fig.4.
20
10
0
y (m)
-10
-20 Anchor
Estimated Trajectory
Test Trajectory
Train Trajectory
-30
Estimated Trajectory
-30 -20 -10 0 10 20 30 40 50
x (m)
Figure 4: The simulation environment and Test1 results
0.6
Test1
0.5
Test2
Test3
0.4 Known Anchor Position
RMSE (m)
0.3
0.2
0.1
0
0 500 1000 1500 2000 2500 3000 3500 4000 4500
5
4
3
GDOP
2
1
0
0 500 1000 1500 2000 2500 3000 3500 4000 4500
trajectory steps
Figure 5: The top figure illustrate RMSE for the four individual testRMSE. The bottom shows GDOP of
test trajectory
5.3. Results of Test Set
In the test set, we can use the anchors estimated in the training set to assist terminal positioning.
Terminal walking through the trajectory shown in Fig.4. We also used the common simulation
parameters described above in our test set. We consider three different parameter settings
dubbed Test1, Test2, and Test3. In Test1 and Test2, we used detection probability pD = 0.9. The
mean number of false alarms are λc = 1 and λc = 2 respectively. In Test3, we used pD = 0.5
and λc = 2 to analyze the robustness of our algorithm to deplorable wireless signal conditions.
We use measurement standard variance σz = 0.01. 10000 particles represented the posterior
pdf of the terminal state. The number of message passing iterations for DA is limited by a
maximum iteration number or the message difference lower than 10−7 .
As a performance benchmark for the accuracy of terminal localization, we also plot in Fig.5
the terminal position RMSE obtained for Test1 with the known anchor real positions. Table ??
shows the 1 − σ cumulative distribution function (CDF) of position error of Test1, Test2, Test3
and known anchors taken together. Test3 is just 0.3m larger than Test2, which suggests a high
accuracy and robustness.
Fig.5 shows terminal position RMSEs of our algorithm obtained individually for the four tests.
We do not give the entire map to plot more detail of four tests. We can find our algorithm plays
a good performance before step 3000. Most errors are less than 0.15. The terminal position error
in Fig.5 illustrates considerable errors between 3500 and 4000. We can find there is a big vibrate
in the bottom left corner. The analysis shows that the anchors that can receive signals in the
bottom right corner are the three vertical anchors, and the other two horizontal anchors cannot
receive signals. We believe that the geometric position of the three vertical anchor points is
responsible for the poor positioning results.
For this reason, we introduce (Geometric Dilution of Precision) GDOP[22] as an analytical
method to analyze measurement errors due to the geometric distribution of anchors. GDOP is
not related to measurement errors but only to the geometric relationship of the anchor, which
fully reflects the geometric affected of anchors’ distribution. The GDOP plot in Fig.5 reflects a
massive error in the lower left-hand corner. At the same time, the fluctuations between 1700
and 2200 also show pool GDOP at the lower right corner. The good thing is that the geometry
of the anchors visible in the bottom right corner is better than that of the anchors in the bottom
left corner, which our algorithm can tolerate. That is why the bottom right corner is not as bad
as the bottom left corner.
6. conclusions and future research direction
In this paper, we proposed a TOA-based anchor mapping and terminal localization algorithm.
The simulation result shows that the FMSC-PHD method based on PHD filtering well estimates
the anchors in the scenarios and solves the problems of DA, false alarm, and miss detection
in indoor localization of wireless signals. In the test set, the factor graph-based BP algorithm
solved both the DA and location estimation problem simultaneously. It was able to be extended
to multi-objective contexts.
Several problems remain to be solved. A limitation of this study is that our search is totally
based on simulation. Therefore, in the future, we are dedicated to validating it with real signals.
Furthermore, our anchor estimation relies on the training set’s observations accuracy. The
test set also depends on the accuracy of the estimation of anchors. We are still working on
PHD methods that can dynamically estimate scenario information, for example, by introducing
information such as AoAs and AoDs for information fusion.
Acknowledgments
This research was funded by National Key Research and Development Program of China with
Grant Number 2018YFB0505200. (Corresponding author: Ao Peng)
References
[1] S. Thrun, Probabilistic robotics, Communications of the ACM 45 (2002) 52–57.
[2] P. S. Farahsari, A. Farahzadi, J. Rezazadeh, A. Bagheri, A survey on indoor positioning
systems for iot-based applications, IEEE Internet of Things Journal (2022).
[3] R. Di Taranto, S. Muppirisetty, R. Raulefs, D. Slock, T. Svensson, H. Wymeersch, Location-
aware communications for 5g networks: How location information can improve scalability,
latency, and robustness of 5g, IEEE Signal Processing Magazine 31 (2014) 102–112.
[4] P. Davidson, R. Piché, A survey of selected indoor positioning methods for smartphones,
IEEE Communications Surveys & Tutorials 19 (2016) 1347–1370.
[5] X. Guo, N. Ansari, F. Hu, Y. Shao, N. R. Elikplim, L. Li, A survey on fusion-based indoor
positioning, IEEE Communications Surveys & Tutorials 22 (2019) 566–594.
[6] Q. D. Vo, P. De, A survey of fingerprint-based outdoor localization, IEEE Communications
Surveys & Tutorials 18 (2015) 491–506.
[7] F. Dwiyasa, M.-H. Lim, A survey of problems and approaches in wireless-based indoor
positioning, in: 2016 International conference on indoor positioning and indoor navigation
(IPIN), IEEE, 2016, pp. 1–7.
[8] W. Liu, Q. Cheng, Z. Deng, H. Chen, X. Fu, X. Zheng, S. Zheng, C. Chen, S. Wang, Survey
on csi-based indoor positioning systems and recent advances, in: 2019 International
Conference on Indoor Positioning and Indoor Navigation (IPIN), IEEE, 2019, pp. 1–8.
[9] C. Gentner, T. Jost, W. Wang, S. Zhang, A. Dammann, U.-C. Fiebig, Multipath assisted
positioning with simultaneous localization and mapping, IEEE Transactions on Wireless
Communications 15 (2016) 6104–6117.
[10] E. Leitinger, P. Meissner, C. Rüdisser, G. Dumphart, K. Witrisal, Evaluation of position-
related information in multipath components for indoor positioning, IEEE Journal on
Selected Areas in communications 33 (2015) 2313–2328.
[11] C. S. Lee, D. E. Clark, J. Salvi, Slam with single cluster phd filters, in: 2012 IEEE International
Conference on Robotics and Automation, IEEE, 2012, pp. 2096–2101.
[12] H. Zhang, S. Y. Tan, Toa based indoor localization and tracking via single-cluster phd
filtering, in: GLOBECOM 2017-2017 IEEE Global Communications Conference, IEEE, 2017,
pp. 1–6.
[13] H. Zhang, S. Y. Tan, C. K. Seow, Toa-based indoor localization and tracking with inaccurate
floor plan map via mrmsc-phd filter, IEEE Sensors Journal 19 (2019) 9869–9882.
[14] R. P. Mahler, Multitarget bayes filtering via first-order multitarget moments, IEEE
Transactions on Aerospace and Electronic systems 39 (2003) 1152–1178.
[15] A. Swain, D. E. Clark, First-moment filters for spatial independent cluster processes, in:
Signal Processing, Sensor Fusion, and Target Recognition XIX, volume 7697, International
Society for Optics and Photonics, 2010, p. 76970I.
[16] A. Swain, D. Clark, The single-group phd filter: An analytic solution, in: 14th International
Conference on Information Fusion, IEEE, 2011, pp. 1–8.
[17] M. Gupta, L. Jin, N. Homma, Static and dynamic neural networks: from fundamentals to
advanced theory, John Wiley & Sons, 2004.
[18] D. Schuhmacher, B.-T. Vo, B.-N. Vo, A consistent metric for performance evaluation of
multi-object filters, IEEE transactions on signal processing 56 (2008) 3447–3457.
[19] J. Williams, R. Lau, Approximate evaluation of marginal association probabilities with
belief propagation, IEEE Transactions on Aerospace and Electronic Systems 50 (2014)
2942–2959.
[20] J. L. Williams, R. A. Lau, Convergence of loopy belief propagation for data association,
in: 2010 Sixth International Conference on Intelligent Sensors, Sensor Networks and
Information Processing, IEEE, 2010, pp. 175–180.
[21] Y. Bar-Shalom, X. R. Li, T. Kirubarajan, Estimation with applications to tracking and
navigation: theory algorithms and software, John Wiley & Sons, 2004.
[22] I. Sharp, K. Yu, Y. J. Guo, Gdop analysis for positioning system design, IEEE Transactions
on Vehicular Technology 58 (2009) 3371–3382.