Open AccessArticle

Analysis and Accuracy Improvement of UWB-TDoA-Based Indoor Positioning System

Paolo Grasso

Mauro S. Innocente

^1,*

Jun Jet Tai

¹,

Olivier Haas

and

Arash M. Dizqah

Autonomous Vehicles & Artificial Intelligence Laboratory (AVAILAB), Centre for Future Transport and Cities, 7th Floor Friars House, Manor House Drive, Coventry CV1 2TE, UK

Centre for Future Transport and Cities, 7th Floor Friars House, Manor House Drive, Coventry CV1 2TE, UK

Smart Vehicles Control Laboratory (SVeCLab), University of Sussex, Brighton BN1 9RH, UK

Author to whom correspondence should be addressed.

Sensors 2022, 22(23), 9136; https://doi.org/10.3390/s22239136

Submission received: 31 October 2022 / Revised: 17 November 2022 / Accepted: 18 November 2022 / Published: 24 November 2022

(This article belongs to the Collection Navigation Systems and Sensors)

Download

Browse Figures

Versions Notes

Abstract

Positioning systems are used in a wide range of applications which require determining the position of an object in space, such as locating and tracking assets, people and goods; assisting navigation systems; and mapping. Indoor Positioning Systems (IPSs) are used where satellite and other outdoor positioning technologies lack precision or fail. Ultra-WideBand (UWB) technology is especially suitable for an IPS, as it operates under high data transfer rates over short distances and at low power densities, although signals tend to be disrupted by various objects. This paper presents a comprehensive study of the precision, failure, and accuracy of 2D IPSs based on UWB technology and a pseudo-range multilateration algorithm using Time Difference of Arrival (TDoA) signals. As a case study, the positioning of a

4 \times 4 m^{2}

area, four anchors (transceivers), and one tag (receiver) are considered using bitcraze’s Loco Positioning System. A Cramér–Rao Lower Bound analysis identifies the convex hull of the anchors as the region with highest precision, taking into account the anisotropic radiation pattern of the anchors’ antennas as opposed to ideal signal distributions, while bifurcation envelopes containing the anchors are defined to bound the regions in which the IPS is predicted to fail. This allows the formulation of a so-called flyable area, defined as the intersection between the convex hull and the region outside the bifurcation envelopes. Finally, the static bias is measured after applying a built-in Extended Kalman Filter (EKF) and mapped using a Radial Basis Function Network (RBFN). A debiasing filter is then developed to improve the accuracy. Findings and developments are experimentally validated, with the IPS observed to fail near the anchors, precision around

\pm 3 cm

, and accuracy improved by about

15 cm

for static and

5 cm

for dynamic measurements, on average.

Keywords:

IPS; ultra-wideband; time difference of arrival; Cramér–Rao lower bound; CRLB; bifurcation curve; debiasing; filtering

1. Introduction

Positioning that is accurate and precise as well as robust and reliable has become an essential part of many applications which require determining the position of an object in space, such as monitoring the location of assets, people, and goods and assisting navigation systems with varying degrees of autonomy while operating within potentially complex and dynamic environments [1,2]

A variety of positioning systems exist which make use of different (i) technologies, (ii) signal properties, and (iii) positioning algorithms. Technologies include inertial navigation systems (INS) [3,4], sound waves [5,6], infrared [7], visible light [8], and radio frequency, including Ultra-Wide Band (UWB) [9], Bluetooth [10], Bluetooth Low Energy (BLE) [11], ZigBee, Wireless Local Area Network (WLAN) [12,13], and Wireless Underground Sensor Network [14]. Signal properties used for positioning include Angle of Arrival (AoA) [7], Time of Arrival (ToA) [14], Time Difference of Arrival (TDoA), Received Signal Strength Indication [15], Time of Flight (ToF), Return Time of Flight (RToF), and Phase of Arrival (PoA) [2]. Positioning algorithms, which include triangulation, trilateration [16,17], proximity, and Two-Way Ranging algorithms [18], are then used in conjunction with the aforementioned technologies and signal properties to estimate and/or track the position of an object. Positioning algorithms tend to be overly sensitive to external disturbances, and therefore often rely on sensor fusion. While this is most often a form of Kalman filter, machine learning models such as neural networks [19,20], clustering algorithms [21], or Bayesian models [22,23] can be used.

Global Navigation Satellite Systems (GNSS) are suitable for efficient outdoor long-range positioning. While the most common technology is the Global Positioning System (GPS), the European Galileo started providing services in 2016 with a constellation of 26 satellites [24]. A GNSS allows an electronic receiver to determine its position by trilateration using radio signal travel times (ToA) from at least four satellites [25]. However, because these signals cannot penetrate walls or objects, use of this technology for Indoor Positioning Systems (IPSs) is infeasible. Conversely, UWB technology is well-suited for IPSs, as they are characterised by large bandwith, high data transfer rates over short distances, short message length, low transmission power, and high obstacle penetration capability [1,9]. In addition, UWB-based IPSs constitute one of the most accurate and precise positioning technologies at present, and are arguably the best choice among current technology [9,26] despite their susceptibility to interference caused by metallic materials or systems working on similar frequencies. Considering the recent drive towards autonomy and self-organisation in robotics (e.g., [27,28]), the precision and accuracy of the IPS is crucial for performing indoor experiments efficiently and safely.

Traditionally, the precision of positioning systems has been studied by performing a Cramér–Rao Lower Bound (CRLB) analysis [29] from the signals perspective, then applying coefficients such as the Geometric Dilution of Precision (GDoP) [30] to capture the geometrical features, e.g., to identify where there is a sudden drop in IPS performance. Examples include the calculation of bounds for ToA [31,32], anchor time synchronisation [33,34,35], and AoA [36]. Such work is of critical importance, as the performance of UWB-based IPSs tends not to be homogenous across the measurement domain [37]; thus, characterising its performance is crucial to understanding the limits of the system and potentially optimising its design and/or layout [38].

CRLB analysis is widely accepted for positioning systems in which the tagged object to be localised/tracked is not in close proximity to the anchors; in such scenarios, localisation algorithms tend to fail due to flipping uncertainty, a well-known problem of geometrical origin [39,40]. Regions in the measurement domain where failure occurs due to geometrical constraints are not considered in the CRLB analysis [38,41]. Furthermore, the possibility of signal degradation due to the anisotropic signal transmission properties of the anchors is almost never considered, even though it is an important property that many antenna optimisation studies have accounted for [42,43,44,45].

This paper is concerned with UWB-based IPSs aimed at localising a tagged moving object (receiver) based on the spatial distribution of the anchors (transceivers) using pseudorange multilateration algorithm and signal TDoA measurements. More specifically, it is concerned with the planar localisation performance within a homogeneous medium with negligible interference and reverberation, while taking into account the anisotropic radiation pattern of an actual DWM1000 antennae module rather than simply assumming ideal signal distributions. A rigorous analysis of system performance is carried out in terms of precision, failure, and accuracy, informed by previous work reported in [38,39,40,46,47,48,49,50]. To summarise, the contributions of this paper are as follows:

(i): A theoretical study of the precision of the position estimates is performed based on a CRLB analysis for round-robin scheduling and an anisotropic representation of the signal-to-noise ratio function of the 3D radiation pattern of the anchor antennas.
(ii): A geometrical study of the 2D IPS domain is carried out, defining bifurcation envelopes that bound the areas where the IPS is predicted to fail. This complements the CRLB analysis, which does not predict regions of failure. Together, they define the so-called flyable area in which positioning is reliable.
(iii): Experiments using an existing IPS with four anchors and a static tagged object are used to validate the precision and failure predictions and to estimate the bias (inaccuracy).
(iv): A debiasing filter is developed to increase the accuracy of the static position estimates, which is then tested for both static and moving tagged objects.

The remainder of this paper is organised as follows: Section 2 presents a comprehensive study of the precision and failure of IPSs based on UWB technology and a pseudorange multilateration algorithm using signal TDoA, including a description of the IPS under study in Section 2.1, a theoretical study of precision using CRLB analysis in Section 2.2, and a geometrical study of failure using a bifurcation envolope analysis in Section 2.3. Section 3 proposes a process to measure the accuracy of the IPS and a debiasing filter to improve it. Section 4 presents the design of the validation and testing experiments, with the results discussed in Section 5. Finally, our conclusions are drawn in Section 6.

2. Theoretical Study of Precision and Failure

The purpose of an IPS is to localise indoor tagged objects (receivers) using the spatial distribution of anchors at known locations (transceivers). The aim here is to study the precision and failure of the system. The former is studied using a CRLB analysis that is specific for TDoA algorithms with round-robin scheduling, while a study of the bifurcation envelope is carried out to identify areas where the system is expected to fail.

2.1. IPS under Study

The system to be studied uses bitcraze’s Loco Positioning System [51] and consists of a drone to be localised and four transceiver anchors positioned at the vertices and facing the centre of a

4 \times 4 m^{2}

domain, as shown in Figure 1. All antennas under study are at a height of 20 cm above the floor; the drone is mounted on a sliding ground-referenced measurement system parallel to the floor equipped with a laser pointer aligned with the onboard UWB antenna to achieve reference positioning with high precision (±1 mm) and accuracy (see Figure 2). This experimental setup is used in Section 4 and Section 5, where the regularly spaced markers on the floor are the sampling positions to be used.

It is important to note that this work assumes that the tagged object to be localised acts strictly as a passive receiver, whereas in the literature it is often a transceiver. Nonetheless, the theoretical results are applicable to both cases as long as the receivers are sensitive and approximately omnidirectional.

2.2. CRLB Analysis for Pseudo-Range Multilateration with Round-Robin Scheduling

The Cramér–Rao Lower Bound (CRLB) analysis is generally deemed suitable for evaluating the precision of an unbiased IPS. It is based on the concept of the Fisher Information Matrix (FIM) for the likelihood of obtaining a correct measurement. For more details on the theory and terminology, refer to Appendix B. The elements of the total FIM for the general positioning problem [48,52] are as shown in Equation (1).

{FIM}_{i j} = {(\frac{\partial h (x)}{\partial x_{i}})}^{T} F_{τ}^{- 1} (x) (\frac{\partial h (x)}{\partial x_{j}}) + \frac{1}{2} t r (F_{τ}^{- 1} (x) \frac{\partial F_{τ} (x)}{\partial x_{i}} F_{τ}^{- 1} (x) \frac{\partial F_{τ} (x)}{\partial x_{j}})

(1)

where

h (x)

is the range vector (e.g., distance between receiver and anchors),

F_{τ}

is the covariance matrix of the

\hat{τ}

measurements, and

t r (\cdot)

is the trace function. Equation (1) considers that the standard deviations (

σ_{i}

) of the likelihood function (and hence

F_{τ}

) vary in space. One column of the Jacobian matrix of

h

is defined as in Equation (2).

\frac{\partial h (x)}{\partial x_{i}} = [\begin{matrix} \frac{\partial h_{12} (x)}{\partial x_{i}} \\ \frac{\partial h_{23} (x)}{\partial x_{i}} \\ \frac{\partial h_{34} (x)}{\partial x_{i}} \\ \frac{\partial h_{41} (x)}{\partial x_{i}} \end{matrix}]

(2)

TDoA measurements for N anchors and round-robin scheduling are referred to as

τ_{r r} = \{τ_{12}, τ_{23}, \dots, τ_{N 1}\}

. Thus, the divergence matrix of

h

for

τ_{r r} = \{τ_{12}, τ_{23}, τ_{34}, τ_{41}\}

for the TDoA2 protocol used by the Loco Positioning System [51,53] is as in Equation (3).

\frac{\partial h (x)}{\partial x} = {[\begin{matrix} \frac{x - x_{1}}{∥x - x_{1}∥} - \frac{x - x_{2}}{∥x - x_{2}∥} & \frac{y - y_{1}}{∥x - x_{1}∥} - \frac{y - y_{2}}{∥x - x_{2}∥} \\ \frac{x - x_{2}}{∥x - x_{2}∥} - \frac{x - x_{3}}{∥x - x_{3}∥} & \frac{y - y_{1}}{∥x - x_{1}∥} - \frac{y - y_{2}}{∥x - x_{2}∥} \\ \frac{x - x_{3}}{∥x - x_{3}∥} - \frac{x - x_{4}}{∥x - x_{4}∥} & \frac{y - y_{3}}{∥x - x_{3}∥} - \frac{y - y_{4}}{∥x - x_{4}∥} \\ \frac{x - x_{4}}{∥x - x_{4}∥} - \frac{x - x_{1}}{∥x - x_{1}∥} & \frac{y - y_{4}}{∥x - x_{4}∥} - \frac{y - y_{1}}{∥x - x_{1}∥} \end{matrix}]}_{4 \times 2}

(3)

Making use of the linear properties of the expected value, the diagonal elements of

F_{τ}

can be calculated as in Equation (4), and its connected elements for consecutive estimators as in Equation (5), where only the

j^{th}

anchor is in common; the hat identifies measurements, the bar stands for the mean, and

E [\cdot]

stands for expectation.

\begin{matrix} F_{i j, i j} = E [{({\hat{τ}}_{i j} - {\bar{τ}}_{i j})}^{2}] = E [{(({\hat{τ}}_{i} - {\bar{τ}}_{i}) - ({\hat{τ}}_{j} - {\bar{τ}}_{j}))}^{2}] = \\ E [{({\hat{τ}}_{i} - {\bar{τ}}_{i})}^{2}] + E [{({\hat{τ}}_{j} - {\bar{τ}}_{j})}^{2}] - 2 E [({\hat{τ}}_{i} - {\bar{τ}}_{i}) ({\hat{τ}}_{j} - {\bar{τ}}_{j})] = σ_{i}^{2} + σ_{j}^{2} \end{matrix}

(4)

\begin{matrix} F_{i j, j k} = E [({\hat{τ}}_{i j} - {\bar{τ}}_{i j}) ({\hat{τ}}_{j k} - {\bar{τ}}_{j k})] = \\ E [(({\hat{τ}}_{i} - {\bar{τ}}_{i}) - ({\hat{τ}}_{j} - {\bar{τ}}_{j})) (({\hat{τ}}_{j} - {\bar{τ}}_{j}) - ({\hat{τ}}_{k} - {\bar{τ}}_{k}))] = \\ E [- {({\hat{τ}}_{j} - {\bar{τ}}_{j})}^{2}] = - E [{({\hat{τ}}_{j} - {\bar{τ}}_{j})}^{2}] = - σ_{j}^{2} \end{matrix}

(5)

Estimating the covariance between seemingly uncorrelated TDoA measurements (

{\hat{τ}}_{i j}

{\hat{τ}}_{k p}

) is not trivial. From the Cauchy–Bunyakovsky–Schwarz inequality, we can derive Equation (6) where

Δ {\hat{τ}}_{i j} = ({\hat{τ}}_{i j} - {\bar{τ}}_{i j})

\begin{matrix} {(E [Δ {\hat{τ}}_{i j} \cdot Δ {\hat{τ}}_{k p}])}^{2} ⩽ E [Δ {\hat{τ}}_{i j}^{2}] \cdot E [Δ {\hat{τ}}_{k p}^{2}] \Rightarrow \\ - \sqrt{E [Δ {\hat{τ}}_{i j}^{2}] \cdot E [Δ {\hat{τ}}_{k p}^{2}]} ⩽ E [Δ {\hat{τ}}_{i j} \cdot Δ {\hat{τ}}_{k p}] ⩽ \sqrt{E [Δ {\hat{τ}}_{i j}^{2}] \cdot E [Δ {\hat{τ}}_{k p}^{2}]} \end{matrix}

(6)

From Equations (5) and (6), the covariance between

{\hat{τ}}_{i j}

and

{\hat{τ}}_{k p}

can be bounded as shown in Equation (7).

\begin{matrix} 0 > F_{i j, k p} = E [({\hat{τ}}_{i j} - {\bar{τ}}_{i j}) ({\hat{τ}}_{k p} - {\bar{τ}}_{k p})] \geq - \sqrt{E [{({\hat{τ}}_{i j} - {\bar{τ}}_{i j})}^{2}] \cdot E [{({\hat{τ}}_{k p} - {\bar{τ}}_{k p})}^{2}]} = \\ - \sqrt{(σ_{i}^{2} + σ_{j}^{2}) (σ_{k}^{2} + σ_{p}^{2})} . \end{matrix}

(7)

Thus, the information matrix of the TDoA measurements set

τ_{r r} = \{τ_{12}, τ_{23}, τ_{34}, τ_{41}\}

for four coplanar anchors using an efficient unbiased estimator is provided by the measurement covariance matrix in Equation (8), where

s_{i}

stands for

σ_{i}^{2}

\begin{matrix} F_{τ} = {[\begin{matrix} s_{1} + s_{2} & - s_{2} & F_{12, 34} & - s_{1} \\ - s_{2} & s_{2} + s_{3} & - s_{2} & F_{23, 41} \\ F_{12, 34} & - s_{3} & s_{3} + s_{4} & - s_{4} \\ - s_{1} & F_{23, 41} & - s_{2} & s_{4} + s_{1} \end{matrix}]}_{4 \times 4} \\ F_{12, 34} = - \sqrt{(s_{1} + s_{2}) (s_{3} + s_{4})} \\ F_{23, 41} = - \sqrt{(s_{2} + s_{3}) (s_{4} + s_{1})} \end{matrix}

(8)

Kaune et al. [48] suggest that the variance for a specific source is as in Equation (9):

\begin{matrix} σ_{i}^{2} (r) = \{\begin{matrix} \frac{a}{{SNR}_{0}} \cdot \frac{r_{i}^{2}}{r_{0}^{2}} & if r_{i} \geq r_{0} \\ \frac{a}{{SNR}_{0}} & if r_{i} < r_{0} \end{matrix} \\ with a = \frac{c^{2}}{B^{2}} \end{matrix}

(9)

where SNR₀ is the signal-to-noise power ratio at a threshold distance

r_{0}

from the

i^{th}

anchor under consideration, c is the signal propagation speed, and B is the bandwidth of the received signal. The SNR₀ varies with the view angle

θ

if the antenna has some directionality. In order to evaluate the SNR(

x

), we use the Friis formula for noise, which provides the relation between the signal gain (over noise) and distance between the transmitter and receiver for different channel frequencies.

2.2.1. Signal-to-Noise Ratio Formulation

The SNR in Equation (10) is the ratio between the power of the signal reaching the receiver (

P_{r}

) and the power of the noise (

P_{N}

SNR (d, θ_{t}, ϕ_{t}, f_{ref}, T, P_{t}) = \frac{P_{r}}{P_{N}}

(10)

It can be expressed as a function of the distance between transmitter and receiver (d), the representative transmission frequency (

f_{ref}

) and bandwidth of the selected channel, the temperature of the environment (T), the transmitting power (

P_{t}

), and the gains of both the transmitting antenna (

G_{t}

) and the receiving antenna (

G_{r}

). Because the receiving antenna is usually very sensitive,

G_{r}

may be neglected in this analysis. The gain

G_{t}

can be a function of the azimuth (

θ_{t}

) and elevation (

ϕ_{t}

) angles with respect to the frame of reference centered on the antenna. The power at the end of the transmission line can be expressed using the contemporary Friis law, as shown in Equation (11):

P_{r} = \frac{P_{t} \cdot G_{t} \cdot G_{r}}{L_{t} \cdot L_{r}} \cdot {(\frac{c}{4 π \cdot f_{ref} \cdot d})}^{2}

(11)

where

L_{t}

and

L_{r}

are the electric losses in the electronics of the transmitter and receiver modules, respectively, which have been embedded in the gains

G_{t}

and

G_{r}

. Here, it is convenient to express everything in logarithmic form.

Combining Equations (9) and (11), the upper bound of the standard deviation is obtained as in Equation (12), where dBm stands for dB milli-watts;

P_{t dBm} (T, V_{i})

is an experimental curve approximating the relationship between the transmission power, the ambient temperature, and the input voltage (

V_{i}

), as in Equation (13) [54]; and

G_{t dBi} (θ_{t}, ϕ_{t})

is the measured transmitting antenna gain with respect to an isotropic antenna, which is a 3D radiation pattern function of the azimuth and elevation angles [55]. Note that the noise power is expanded into a thermal noise power term,

k_{B} T B_{w}

, where

k_{B}

is the Boltzmann constant for radiation of a black body (≈

1.38 \times 10^{- 23} J / K

\{\begin{matrix} \begin{matrix} {SNR}_{dB} = P_{t dBm} (T, V_{i}) + G_{t dBi} (θ_{t}, ϕ_{t}) - 10 {log}_{10} (k_{B} T B_{w} 10^{3}) - 20 {log}_{10} (4 π f_{ref} d / c) \end{matrix} \\ σ^{2} = \frac{c^{2}}{B_{w}^{2}} \cdot 10^{- {SNR}_{dB} / 10} \end{matrix}

(12)

P_{t dBm} (T, V_{i}) = P_{t 0 dBm} + {\frac{\partial P_{t}}{\partial T}|}_{T_{ref}} (T - T_{ref}) + {\frac{\partial P_{t}}{\partial V}|}_{V_{ref}} (V_{i} - V_{ref})

(13)

2.2.2. Radiation Pattern of the DW1000 Anchor Antenna

Most theoretical studies tend to use ideal distributions of the signal around the antennae (isotropic, bi-conical, cardioid, unidirectional, etc.); however, we hypothesise that the actual signal distribution is important, as it may have a non-negligible impact on the precision of the IPS.

In order to reconstruct the 3D radiation pattern from the three measured sections in the azimuth (

θ)

, elevation-1 (

ϕ_{1}

) and elevation-2 (

ϕ_{2}

) planes (see Figure 3), we formulate a linear combination of the boundary values of the considered quadrant. Using the system of Equations (14), the 3D radiation pattern depicted in Figure 4 can be obtained.

\{\begin{matrix} a_{1} = c o s^{2} (θ) \cdot (1 - c o s^{40} (ϕ)) \\ a_{2} = (1 - c o s^{2} (θ)) \cdot (1 - c o s^{40} (ϕ)) \\ a_{3} = c o s^{40} (ϕ) \\ G (θ, ϕ) = a_{1} \cdot G_{ϕ_{1}} + a_{2} \cdot G_{ϕ_{2}} + a_{3} \cdot G_{θ} \end{matrix}

(14)

The obtained radiation pattern of the antennas supports our choice to place the anchors in our experiments facing the centre of the domain (see Section 2.1), despite BitCraze seemingly recommending that they be placed in pairs facing each other and forming a 90-degree angle with the floor (see figure with eight anchors placed in a room in [56]). Furthermore, this study provides additional variables for the optimal design problem formulation of IPSs, namely, the antenna orientations.

2.2.3. Analytical Results of CRLB Analysis

The CRLB analysis carried out here considers two different representative distributions of four anchors, namely, a symmetrical layout and a random one, as shown in Figure 5. The ripples of the contour lines in Figure 5 are to be expected due to the anisotropy of the radiation pattern in Figure 4.

The best precision is obtained within the convex hull of the anchors. In this case, it is about

\pm 5

cm with 99% confidence level (i.e.,

k = 2.58

). A realistic non-isotropic transmitting antenna gain (DWM1000 module [55]) is applied for the estimation of the SNR, hence the slight fluctuations in the represented values.

2.3. Bifurcation Envelope Analysis

The CRLB analysis is crucial for estimating the precision of the positioning system across the domain of measurements, although the analysis does not consider regions in the measurement domain where failure occurs due to geometrical constraints [38,41]. For instance, IPSs suffer from a well-known problem of geometrical origin called flipping uncertainty [39,40]. A TDoA map, which is a geometrical representation of the TDoA measurements, can be used to address this issue. Thus, we define a so-called flyable area using a combination of the CRLB analysis carried out in Section 2.2 and the bifuraction envelope derived from a geometrical study carried out in this section. Specifically, this flyable area defines a region of space within which system precision is guaranteed to be inside the bounds calculated using the CRLB analysis. It follows that any object’s location measured should only be trusted when within this domain.

2.3.1. Bifurcation Curve

The bifurcation curve is the projection of the TDoA map boundaries from the

τ

-plane (pseudo-range space) to the space of source localisation (2D in this case). The bifurcation curve, as defined in [39], is the quintic curve

\tilde{E} (x)

depicted by the roots of a polynomial

P (x)

which represents the TDoA map constraints. The definition of

P (x)

and examples of algebraic equations of

\tilde{E} (x)

can be found in [47], while its rigorous derivation is presented in [39] using tools such as exterior algebra formalism and Minkowski space. This formulation is invariant under permutation of the TDoA measurements, which means that scheduling does not affect this analysis. Any TDoA-based system has a unique solution of the positioning problem if

P (x)

is negative, which defines the region outside the bifurcation curves surrounding the anchors. The multilateration algorithm within the bifurcation curves (convex regions) returns either two mirrored solutions or complex solutions with no physical meaning. An example of a bifurcation curve is shown in Figure 6a,b for the case of three anchors {

m_{2}, m_{3}, m_{4}

2.3.2. Bifurcation Envelope

For positioning systems comprising several antennas, the bifurcation curve changes dynamically depending on the paired times of arrival (TOAs) considered in each TDoA query. As discussed earlier, the system fails to estimate the position of an object within the concave regions of the bifurcation curves (i.e., those containing the anchors). In order to ensure a unique solution for any possible pairing, a so-called bifurcation envelope is defined which bounds all bifurcation curves on each anchor (e.g., one curve surrounding each anchor for three anchors and four curves for four anchors). In Figure 6c,d, the flyable area shaded in yellow is defined as the intersection of two areas:

1.: The unique-solution area, defined as the intersection of all concave areas outside each green bifurcation envelope (i.e., not including anchors).
2.: The region with acceptable precision returned by the CRLB analysis (the convex hull).

The

i^{th}

transmitting anchor is represented by

m_{i}

, with

m_{1}

being disregarded in Figure 6a. The centroids of triplets (

m_{i}

m_{j}

m_{k}

) are represented by

C_{i j k}

, while C is the collective centroid. Figure 6b,c shows the flyable area (shaded in yellow) and the convex hull of the four anchors (dotted magenta trapezoid, acceptable precision).

3. Filter Design

Thus far, we have analysed the precision and failure of UWB-based IPSs based on a pseudo-range multilateration algorithm with round-robin scheduling and signal TDoA. The aim in this section is to develop a filtering process to improve the accuracy of the system.

3.1. Proposed Filter Design

In our experimental setting, the object to be positioned is a Crazyflie 2.0 nano-quadcopter and the IPS is bitcraze’s Loco Positioning System [56]. This is setup already equipped with an Extended Kalman Filter (EKF) [57,58], which transforms raw sensor measurements into better estimates of the state of the drone (i.e., higher precision). The EKF developers note that the position estimates are affected by a measurement bias which appears to be non-uniform in space. In other words, the quadcopter is estimated to be in a position that is shifted from the actual one. To address this issue, we proposed that a debiasing filter be incorporated here after the built-in EKF.

In addition to a plain EKF, other practitioners may wish to incorporate additional filters to improve precision. An attempt to do this is discussed in Appendix C, although no meaningful precision improvement was observed with that particular set of filters.

3.2. Debiasing Filter

In this section, a filter is proposed and developed aiming to reduce systematic biases. Specifically, the aim of this Debiasing Filter (DF) is to increase the accuracy of the localisation of the drone by subtracting the expected bias from the measurements. Assuming that discrete distributions of variances and biases (see Section 4.1) have been obtained by statistical post-processing of consecutive position measurements, two major problems arise:

1.: The bias values are available only at a limited set of points, and therefore they need to be interpolated to cover the continuous domain.
2.: The bias to be subtracted from a measured position to obtain the actual one is a function of the actual position itself.

To address the first problem, a continuous model must be fitted to the limited data. This takes the form of a surface for 2D positioning (not necessarily defined on a regular grid) and a hypersurface in 4D space for 3D positioning. To address the second problem, a bias map must be built as a function of the measured rather than the actual positions.

In order to explain the proposed debiasing process, we begin by defining the debiased measurement as

\hat{x}

, the measured posistion as

\bar{x}

, and the cloud of actual positions as

X_{i j}

, as shown in Equation (15) and Figure 7.

\begin{array}{l} \begin{matrix} x = (x, y) : \end{matrix} & P o s i t i o n \\ \bar{x} : & M e a s u r e d x \\ X_{i j} = (X_{i j}, Y_{i j}) : & C l o u d o f a c t u a l x c o r r e s p o n d i n g t o (i, j) p o i n t \\ \hat{x} : & D e b i a s e d x \end{array}

(15)

Any measured position can be expressed as the sum of the actual position, a bias value (b), and a noise value (

R

), as shown in Equation (16) for independent x and y components. Note that both b and

R

depend on the actual position, where p stands for the variance.

\begin{matrix} \bar{x} = X_{i j} +^{x} b (X_{i j}) + R (^{x} p (X_{i j}), x) \\ \bar{y} = Y_{i j} +^{y} b (X_{i j}) + R (^{y} p (X_{i j}), x) \end{matrix}

(16)

Assuming

R

to be negligible, the real position can be obtained from the measurement

\bar{x}

by simply subtracting the bias associated with the real position itself (see Equation (16) and Figure 7). It would be more practical to have the bias as a function of the measured positions rather than the actual positions. To this end, a set of bias measurements is obtained at a regular grid of points of known locations and added to the positions where they were measured, thereby obtaining an irregular grid. By fitting a model to the bias measurements associated with their corresponding points in the new grid, a bias map is obtained as a function of the measured rather than the actual postions:

^{x} β (\bar{x})

and

^{y} β (\bar{x})

For the purpose of the following derivation, continuous interpolating functions of the x and y biases, both along the

\hat{i}

and

\hat{j}

axes, have to be obtained (e.g., cubic splines). Therefore, from the biases of x measurements around any

x_{i j}

position, two interpolating functions can be obtained, namely,

_{\hat{i}}^{x} b_{i j}

and

_{\hat{j}}^{x} b_{i j}

along the

\hat{i}

and

\hat{j}

axes, respectively. Likewise,

_{\hat{i}}^{y} b_{i j}

and

_{\hat{j}}^{y} b_{i j}

interpolating functions can be defined for the biases of y measurements. The aim is to write weighted averages of the biases around the estimation position in order to estimate the expected biases. However, instead of performing a surface integral, the average of two integrals in perpendicular directions is considered. Figure 8 (left) depicts the problem in the

\hat{i}

direction for the bias of the measurement of the x-component of the position. Thus, the interpolated bias function

_{\hat{i}}^{x} b_{i j}

multiplied by the weighting probability distribution

^{x} γ_{i j}

is integrated in the x direction and normalised by the length of the considered interval. A coverage factor of

k = 3

is set, which means that approximately 99% (level of confidence) of the measurements of the real position

{^{r} x}_{i j}

fall in the interval between

{(x - 3^{x} σ)}_{i j}

and

{(x + 3^{x} σ)}_{i j}

, where

^{x} σ_{i j}

is the standard deviation of the Gaussian distribution

^{x} γ_{i j}

. The same integral can be evaluated in the

\hat{j}

direction and the two integral values can be averaged in order to obtain the corrected bias value of the x-component measurements, as shown in Equation (17) and Figure 8 (right). The same process can be applied for evaluating the corrected bias of the y-component measurements, as in Equation (18).

^{x} b_{i j}^{*} = \frac{1}{2} \int_{- 3 \cdot^{x} σ_{i j}}^{+ 3 \cdot^{x} σ_{i j}}^{x} γ_{i j} \cdot_{\hat{i}}^{x} b_{i j} \cdot d x + \frac{1}{2} \int_{- 3 \cdot^{y} σ_{i j}}^{+ 3 \cdot^{y} σ_{i j}}^{y} γ_{i j} \cdot_{\hat{j}}^{x} b_{i j} \cdot d y

(17)

^{y} b_{i j}^{*} = \frac{1}{2} \int_{- 3 \cdot^{x} σ_{i j}}^{+ 3 \cdot^{x} σ_{i j}}^{x} γ_{i j} \cdot_{\hat{i}}^{y} b_{i j} \cdot d x + \frac{1}{2} \int_{- 3 \cdot^{y} σ_{i j}}^{+ 3 \cdot^{y} σ_{i j}}^{y} γ_{i j} \cdot_{\hat{j}}^{y} b_{i j} \cdot d y

(18)

For the remaining derivations, refer to Figure 9. By rewriting the decomposition shown in Figure 7 and neglecting the noise component, the measured position is shifted from the original one approximately by the corrected weighted biases expressed in Equations (17) and (18) as shown in Equation (19).

\begin{matrix} {\bar{X}}_{i j} = X_{i j} +^{x} b_{i j}^{*} \\ {\bar{Y}}_{i j} = Y_{i j} +^{y} b_{i j}^{*} \end{matrix}

(19)

Therefore, while the original experimentally-obtained biases were distributed on a regular quadrangular grid

[{^{r} x}_{i, j} {^{r} y}_{i, j}]

, the new corrected biases

b^{*}

can be distributed over a deformed grid

[{\bar{x}}_{i, j} {\bar{y}}_{i, j}]

. As shown in Equation (20), the two new corrected bias distributions of the x (

^{x} m_{i, j}

) and y (

^{y} m_{i, j}

) measurements can be interpolated, obtaining bias surfaces that are functions of the measured positions.

\begin{matrix} ^{x} m_{i j} = {[\begin{matrix} {\bar{X}}_{i j} & {\bar{Y}}_{i j} & ^{x} b_{i j}^{*} \end{matrix}]}^{T} & \overset{i n t e r p .}{\to}^{x} β (x) \\ ^{y} m_{i j} = {[\begin{matrix} {\bar{X}}_{i j} & {\bar{Y}}_{i j} & ^{y} b_{i j}^{*} \end{matrix}]}^{T} & \overset{i n t e r p .}{\to}^{y} β (x) \end{matrix}

(20)

Finally, it is possible to subtract these new interpolating bias functions from the measured position in order to obtain a debiased measurement, as in Equation (21).

\begin{matrix} \hat{x} = x -^{x} β (x) \\ \hat{y} = y -^{y} β (x) \end{matrix}

(21)

Two examples of calibrated debiasing function fields can be found in Section 4.2 while the formulation of the Radial Basis Function Network (RBFN) used for interpolation of the debiasing values is explained in Appendix A.

4. Design of Experiments

4.1. IPS Bias Map Generation

For simplicity, our version of the IPS (with built-in EKF + DF) is referred to as IPS-2, while the original version (only with built-in EKF) is referred to as IPS-1. In order to build the maps, a large number of measurements (

N = 700

) are taken at a sampling frequency of 100 Hz while keeping the drone stationary for at least 30 s on each marker (

X_{i j}

). The drone is kept aligned with the x axis and parallel to the floor, as the effect of its attitude is not being investigated. Finally, the bias (b), standard deviation (

σ

), and mean squared error (MSE) are computed. For instance, their values in the x direction (superscript

^{x}

) are as in Equation (22), where

x^{(k)}

is the

k^{th}

position measurement. The resulting maps can be found in Figure 10 and Figure 11.

\begin{matrix} ^{x} b_{i j} = N^{- 1} \sum_{k = 1}^{N} x^{(k)} - X_{i j} \\ ^{x} σ_{i j} = {({^{x} MSE}_{i j} - {^{x} b}_{i j}^{2})}^{0.5} \\ ^{x} {MSE}_{i j} = N^{- 1} \sum_{k = 1}^{N} {(x^{(k)} - X_{i j})}^{2} \end{matrix}

(22)

4.2. DF Calibration and Validation Setup

Using the obtained bias measurements, the debiasing filter for IPS-2 is calibrated using Radial Basis Function Networks (RBFN); refer to Appendix A for the formulation. This produces a final output similar to Figure 12 or Figure A2.

For validation, the estimates provided by IPS-2 are compared to those returned by IPS-1 at a predefined set of points (shown in Figure 1) not been previously used for calibration purposes. The variance and bias are evaluated to provide statistical insight into the performance of the newly developed filter.

4.3. DF Validation under Dynamic Setup Conditions

In addition to static validation, the estimates provided by IPS-1 and IPS-2 are further compared in a dynamic system with the vehicle cruising at different speeds.

The drone is mounted on a mobile stand which is constrained to move along an encoder rail (see Figure 2). In the frame of reference along the rail, the position (s) of the drone, and thus that of the optical sensor, is estimated as in Equation (23):

s^{(k)} = Δ s \cdot n^{(k)}

(23)

where

Δ s

is the constant distance between consecutive pins (optical obstacles), n is the counted number of pins, and k is the measurement frequency at which the data from the sensors are recorded. Note that n is supposed to increase more slowly than k.

Figure 13 represents rail and IPS position measurements taken on a horizontal rail. The real position of the drone along the rail (

x_{r}^{(k)}

) in the inertial frame of reference can be obtained by projecting

s^{(k)}

along the angle between the rail and the x-axis,

θ_{r}

\begin{matrix} x_{r}^{(k)} = (x_{r}^{(k)}, y_{r}^{(k)}) \\ x_{r}^{(k)} = s^{(k)} \cdot cos (θ_{r}) and y_{r}^{(k)} = s^{(k)} \cdot sin (θ_{r}) \\ θ_{r} = atan (\frac{y_{b} - y_{a}}{x_{b} - x_{a}}) \end{matrix}

(24)

Every time a new pin is detected at time step p, the estimated position of both IPS-1 and IPS-2 (

x_{1}^{(k)}

x_{2}^{(k)}

) and the actual position on the rail are recorded as variables

ξ

and

υ

\begin{matrix} if n^{(k)} > n^{(k - 1)} then : \\ ξ_{r}^{(p)} = x_{r}^{(k)} υ_{r}^{(p)} = y_{r}^{(k)} τ^{(p)} = t^{(k)} \\ ξ_{1}^{(p)} = x_{1}^{(k)} υ_{1}^{(p)} = y_{1}^{(k)} ξ_{2}^{(p)} = x_{2}^{(k)} υ_{2}^{(p)} = y_{2}^{(k)} \end{matrix}

(25)

In addition, the velocity estimations with IPS-1 and IPS-2 are compared to the discrete average velocity on the rail:

\begin{matrix} ν_{x}^{(p)} = \frac{Δ ξ_{r}^{(p)}}{Δ τ^{(p)}} and ν_{y}^{(p)} = \frac{Δ υ_{r}^{(p)}}{Δ τ^{(p)}} \\ Δ ξ_{r}^{(p)} = Δ s \cdot cos (θ_{r}) = const . \\ Δ υ_{r}^{(p)} = Δ s \cdot sin (θ_{r}) = const . \\ Δ τ^{(p)} = τ^{(p)} - τ^{(p - 1)} \end{matrix}

(26)

where

Δ τ^{(p)}

is the time passed between the detection of the

p^{th}

and

{(p - 1)}^{th}

pins.

Because the aim is to evaluate how well the IPS-2 performs with respect to IPS-1 at different crusing speeds of the drone, multiple measurements are required at different speeds. The recorded positional data are classified in different groups.

4.4. Square Path Experiment Setup

The aim in this experiment is to partially reproduce the experiment with a drone following a square path [57] (Figure 14a). Because we want to investigate only the performance of the debiasing filter of IPS-2, we disable drone flight in favour of it being driven around by the mobile support along the square path. This decouples the performance of the debiasing filter from the control loop of the drone. The estimated positions of IPS-1 and IPS-2 are then compared. The expected result is depicted in Figure 14b.

5. Results and Discussion

5.1. Proof of Accuracy Improvement

The results of the experiments in Section 4.2 are presented in this section along with the accompanying Figure 15 and Figure 16. In order to highlight the overall accuracy gain using the DF, the absolute value of the bias is represented. Note that the calibration points on the main grid are spaced 50 cm apart, while the validation mesh is staggered by 25 cm from the main calibration points.

The whiter areas in Figure 15 and Figure 16 correspond to more accurate areas. It can be seen that the contribution of the proposed DF is evident. It is noticeable that the DF fails to improve the accuracy in a few validation points, which means that the sampling points used for the mapping did not capture the gradient of the bias. A finer sampling mesh would most likely solve this issue. However, a compromise must be made between mapping refinement and the complexity of RBFN interpolation.

Another interesting aspect of the theoretical analysis previously carried out is manifested around the anchor positioned at (0,0). The points within the 50 cm radius around this anchor are undefined because these locations fall within the bifurcation envelope (refer to Section 2.3). No position can be measured in this area, and the DF is expected to fail.

5.2. Dynamic Validation of Debiasing

In this section, the results of the dynamic experiments in Section 4.3 and Section 4.4 are discussed. A reduction in the performance of the DF is expected, partly due to the intrinsic effects introduced by the positioning algorithm that are not addressed by the DF. Though less impressively than for the static case, the proposed DF improves the accuracy of the position estimates, as can be observed in Table 1 and Table 2. More precisely, Table 1 refers to the dynamic validation of the DF explained in Section 4.3, while Table 2 shows the statistics of the results of the square path experiment (Section 4.4).

An example of results from a single dynamic experiment are shown in Figure 17. This experiment was repeated ten times for each rail position in order to obtain the general trend of the IPS measurements.

In Figure 17, note the dynamic misbehaviour of the IPS at around 22 s, which cannot be addressed by the proposed DF, although it could perhaps be handled by other filters (refer to Appendix C). Figure 18 presents a graphical representation of the overall results of the square-path experiment, while Table 2 presents a quantitative summary.

The Root Mean Square Errors (RMSEs) in Table 2 are computed by comparing the trendlines for each edge to the actual position of the drone on the rail at every time step. For instance, the RMSE_IPS-1 on the bottom edge is obtained using the data cloud (cyan colour in Figure 18) of ten experiments on the edge from

(0.5, 0.5)

(3.5, 0.5)

. The same dynamic issues that were pointed out in the validation experiment in Figure 17 persist in the square-path experiment (yellow regions in Figure 18). Assuming that the DF is insufficient to address the intrinsic problems of the IPS under study, we believe it may be useful to isolate this misbehaviour and provide statistics on the data unaffected by this (hence, the raw and sel. columns in Table 2).

6. Conclusions and Future Work

While the precision of IPSs is generally well studied, reasonably estimated, and provided by the manufacturer, failure and accuracy tend to be assessed poorly, if not plainly disregarded. In this paper, we carried out a comprehensive study of the precision, failure, and accuracy of 2D IPSs based on UWB technology and pseudo-range multilateration algorithm with round-robin scheduling using signal TDoA, in addition to developing a debiasing filter to improve accuracy. Although a number of aspects of this investigation are either general or can be generalised, the focus was on a specific setting of the IPS.

A theoretical study of the precision of the position estimates was performed based on a CRLB analysis taking into account the anisotropic radiation pattern of the anchors antennas. The precision is found to be higher within the convex hull of the anchors. Furthermore, visual inspection of the radiation pattern indicates that the orientation of the anchor antennas should probably be towards the centre of the domain.

A geometrical study of the two-dimensional positioning domain was carried out by bounding the areas in which the IPS is predicted to fail via bifurcation envelopes. The intersection between the convex hull and the region outside the bifurcation envelopes results in what we call the flyable area, within which the performance of the IPS is reliable.

The accuracy of the system was measured experimentally on a regular grid of points of known locations after applying a built-in EKF, thereby building a so-called bias map by fitting a Radial Basis Function Network (RBFN). In order to improve accuracy, a debiasing filter was developed by correcting the bias map to ensure that it depends on the measured rather than the actual positions. In this way, it can simply be subtracted from the measurements to debias them.

Our findings and developments were experimentally validated, with the IPS observed to fail near the anchors, precision found to be about

\pm 3 cm

, and accuracy improved by about

15 cm

for static and by

5 cm

for dynamic measurements on average. The proposed method to define the flyable area and build the precision maps, accuracy maps, and debiasing filter is generalisable and repeatable for any other IPS that uses UWB technology and a multilateration algorithm based on the TDoA signal property. The numerical values reported here correspond to the specific IPS used in the experiments.

Future work might involve the generalisation of this study to 3D IPSs, the automation of the process to make it more easily applied in different settings, and the optimisation of the number, positioning, and orientation of the anchors in 3D to account for the slight anisotropy of the radiation pattern of the UWB module.

Author Contributions

The contributions of the authors to the presented research are as follows: conceptualization, P.G., M.S.I. and A.M.D.; methodology, P.G. and M.S.I.; software, P.G.; validation, P.G. and M.S.I.; formal analysis, P.G., M.S.I. and J.J.T.; investigation, P.G.; resources, M.S.I.; data curation, P.G.; writing—original draft preparation, P.G. and M.S.I.; writing—review and editing, P.G., M.S.I., J.J.T., O.H. and A.M.D.; visualization, P.G.; supervision, M.S.I. and O.H.; project administration, M.S.I. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

IPS	Indoor Positioning System
ToA	Time of Arrival
TDoA	Time Difference of Arrival
UWB	Ultra-Wideband
CRLB	Cramér–Rao Lower Bound
GNSS	Global Navigation Satellite System
GDoP	Geometric Dilution of Precision
CRLB	Cramér–Rao Lower Bound
EKF	Extended Kalman Filter
DF	Debiasing filter

Appendix A. Radial Basis Function Network Implementation

In this section, the formulation of a Radial Basis Function Network (RBFN) is introduced for the interpolation of the debiasing distributions in both 2D environments (see Equation (20)). In the following equations, the change of variable expressed in Equation (A1) must be considered to ensure that the RBFN works on strictly positive values.

\begin{matrix} b_{i j} =^{x} b_{i j}^{*},^{y} b_{i j}^{*}, or^{z} b_{i j}^{*} and \\ if b_{\min} = \min \{b_{i j}\} < 0 then b_{i j} = b_{i j}^{*} - b_{\min} \end{matrix}

(A1)

The Gaussian activation function (

f_{i j}

) is defined on every marker (

X_{i j}

) that was previously used for mapping purposes. Hence, every marker represents a node in the network.

\begin{matrix} f_{i j} (x) = c_{1} \cdot b_{i j} \cdot e^{- c_{2} \cdot d^{2}} \\ d^{2} = {∥ x - X_{i j} ∥}^{2} \end{matrix}

(A2)

While the constant

c_{2}

can be uniquely defined for each node (see Equations (A3)–(A6)), the constant

c_{1}

is used for calibration of the RBFN interpolation. The value of

c_{1}

can be chosen between 0 and 1. By trial and error, the best acceptable value was found to be

c_{1} = 0.5

\begin{matrix} f_{i j} (X_{i j}) = c_{1} \cdot b_{i j} \\ f_{i j} (X_{m n}) = c_{1} \cdot b_{i j} \cdot e^{- c_{2} \cdot d_{m n}^{2}} = \frac{1 - c_{1}}{8} b_{m n} \\ d_{m n}^{2} = {∥ X_{m n} - X_{i j} ∥}^{2} \end{matrix}

(A3)

\begin{matrix} ln (c_{3} \frac{b_{m n}}{b_{i j}}) = - c_{2} \cdot d_{m n}^{2} \\ c_{3} = \frac{1 - c_{1}}{8 \cdot c_{1}} \end{matrix}

(A4)

c_{2, m n} = - \frac{1}{d_{m n}^{2}} ln (c_{3} \frac{b_{m n}}{b_{i j}})

(A5)

The constant

c_{2, i j} (θ)

is a periodic function interpolating the values of

c_{2, m n}

on a neighbourhood of eight surrounding points (in 2D).

\begin{matrix} c_{2, i j} (θ) = interp (c_{2, m n}, θ_{m n}) periodic in [- π, + π] \\ m n = [(i + 1, j) (i + 1, j + 1) (i, j + 1) (i - 1, j + 1) (i - 1, j) \\ (i - 1, j - 1) (i, j - 1) (i + 1, j - 1) (i + 1, j)] \\ with θ_{m n} = atan (\frac{Y m n - Y i j}{X_{m n} - X_{i j}}) \end{matrix}

(A6)

Therefore, after all the

c_{2, i j} (θ)

constants have been defined, the resulting RBFN (shown in Figure A1 and Figure A2) is as provided in Equation (A7).

\begin{matrix} b^{(NN)} (x) = b_{\min} + \sum_{i, j}^{N, M} f_{i j} (x) \\ f_{i j} (x) = c_{1} \cdot b_{i j} \cdot e^{- c_{2, i j} \cdot d_{i j}^{2}} \\ d_{i j}^{2} = {∥ x - X_{i j} ∥}^{2} \\ c_{1} = 0.5 \\ c_{2, i j} = c_{2, i j} (θ) periodic in [- π, + π] \\ view angle θ_{i j} = atan (\frac{y - Y i j}{x - X_{i j}}) \end{matrix}

(A7)

Figure A1. RBFN of bias values on marker points (red stars) for estimating the position of the x (left) and y (right) components.

Because the Gaussian activation functions are defined here such that they have influence only on the adjacent nodes, the debiasing filter (DF) code uploaded in the Crazyflie firmware takes account of this in order to optimise the performance in terms of real-time computation. Therefore, not all

f_{i j} (x)

are computed at every time step, only those which are within a distance (

d_{i j}

) of 1 m from the measured position.

b^{(NN *)} (x) = b_{\min} + \overset{N, M}{\sum_{i, j ∣ d_{i j} < 1 m}} f_{i j} (x)

(A8)

Figure A2. RBFN surface interpolating bias values of estimated position of the x (left) and y (right) components.

Appendix B. Reference CRLB Analysis

Before proceeding with the analysis, the reference theory is used to set the background nomenclature. Considering an indoor positioning system consisting of only three anchors and using Time Difference of Arrival (TDoA) multilateration, the pseudo-ranges can be defined as the range differences between the node and each anchor:

\begin{matrix} pseudorange : & τ_{i j} (x) = d_{j} (x) - d_{i} (x) \\ range : & d_{i} (x) = ∥x - x_{i}∥ and i, j = 1, 2, 3 \end{matrix}

(A9)

where

d_{i}

is the distance between the drone (

x

) and the

i^{th}

anchor position (

x_{i}

). According to [39], without loss of generality, the speed of propagation in the medium is considered to be 1 because the range is linearly dependent on the Time of Arrival (ToA). The value of

τ_{i j}

can be obtained by defining a TDoA mapping that transforms the two-dimensional space of the source location to a space of pseudo-ranges called the

τ

-plane, as suggested in [46]:

\begin{matrix} τ_{2} : & R^{2} & \to & R^{2} \\ x & \to & (τ_{12} (x), τ_{13} (x)) \end{matrix}

(A10)

Studying the TDoA map is crucial for the mathematical characterisation of the localisation problem. Considering an IPS that consists of a network of N anchors and one node (e.g., a drone), M measurements (number depending on the localisation algorithm) are obtained at every time step according to the frequency at which messages are sent from the anchors to the node. Every measurement is modeled as a normal distribution which is a function of both the real measurements and additive Gaussian noise. The standard deviation of the noise changes in space with the distance from any transmitting anchor. The collection of M measured pseudo-ranges

\hat{τ}

at time step k can be expressed as

\begin{matrix} {\hat{τ}}^{(k)} \in R^{M}, \\ {\hat{τ}}_{i j} \sim N (τ_{i j} (x), {\bar{σ}}_{i j}^{2}), {\bar{σ}}_{i j} = f (σ_{i}, σ_{j}), \\ τ_{i j} = ∥x - x_{i}∥ - ∥x - x_{j}∥, \\ i, j \in \{1, \dots, N\} with i \neq j, \end{matrix}

(A11)

where

{\hat{τ}}_{i j}

is the individual pseudo-range measurement using the two ToAs of the node to the ith and jth anchors, and

τ_{i j}

is the real range difference. While the superscript

^{(k)}

is omitted, it is implicitly inferred in further analyses. Please note that

\hat{τ}

is a column vector, not a matrix. Here,

{\hat{τ}}_{i j}

can have as many components as the binomial coefficient

(\binom{N}{3}) = \frac{N!}{3! (N - 3)!}

. The functions evaluating the combined standard deviations are presented for the specific IPS under study in Section 2.

The well-known Cramér–Rao Lower Bound (CRLB) analysis is very useful for evaluating the precision of an unbiased IPS. Such analysis is based on the concept of the Fisher Information Matrix (FIM), which contains the likelihood of obtaining a correct measurement. The elements of the total FIM for the general positioning problem according to [48,52] are:

\begin{matrix} {FIM}_{i j} = {(\frac{\partial τ (x)}{\partial x_{i}})}^{T} F_{τ}^{- 1} (x) (\frac{\partial τ (x)}{\partial x_{j}}) \\ + \frac{1}{2} t r (F_{τ}^{- 1} (x) \frac{\partial F_{τ} (x)}{\partial x_{i}} F_{τ}^{- 1} (x) \frac{\partial F_{τ} (x)}{\partial x_{j}}) \end{matrix}

(A12)

where x is the index of

\hat{τ}

measurements and

t r (M)

is the trace of a matrix M. Because each standard deviation is considered to be changing in space, i.e.,

σ_{i} (x, y)

, the correction term (the second row in Equation (A12)) is acknowledged in the following analysis. The likelihood function,

L

, which describes the relative odds of obtaining the observed data

h

for all permissible values of the parameter

x

for a single measurement h, is:

L (\hat{h} | x) = \frac{1}{\sqrt{2 π} σ (x)} e^{(- \frac{1}{2 σ^{2} (x)} {(\hat{h} - h (x))}^{2})}

(A13)

In our positioning problem,

h (x)

is the range, or similarly the ToA, and the parameter

x

is the node position. The standard deviation of this distribution is again

σ

. Considering the Fisher Information Matrix (FIM) in terms of likelihood [48,52], the logarithmic likelihood has to be considered:

ℓ (\hat{h} | x) = ln (\frac{1}{\sqrt{2 π} σ (x)}) - \frac{1}{2} \frac{{(\hat{h} - h (x))}^{2}}{σ^{2} (x)}

(A14)

Hence, the total FIM for the location, considering that each standard deviation changes in space (

σ_{i} = ξ (x, y)

), is as follows:

{FIM}_{i j} = {(\frac{\partial h (x)}{\partial x_{i}})}^{T} F_{τ}^{- 1} (x) (\frac{\partial h (x)}{\partial x_{j}}) + \frac{1}{2} t r (F_{τ}^{- 1} (x) \frac{\partial F_{τ} (x)}{\partial x_{i}} F_{τ}^{- 1} (x) \frac{\partial F_{τ} (x)}{\partial x_{j}})

(A15)

where

t r (M)

is the trace of the matrix M.

Here,

F_{τ}

is the information matrix of the selected set of TDoA measurements. Suppose that the TDoA protocol requires M measurements; then,

F_{τ}

F_{τ} = {[F_{i j}]}_{M \times M}

(A16)

Using an efficient unbiased estimator, it is proven [52] that

F_{τ}

is the measurement covariance matrix:

F_{τ} = {[E [({\hat{τ}}_{i j} - {\bar{τ}}_{i j}) ({\hat{τ}}_{k p} - {\bar{τ}}_{k p})]]}_{M \times M}

(A17)

where the indices

i, j, k, p

depend on the selected scheduling.

As an example, if all the TDoA measurements are performed while keeping anchor 1 as the reference (as in all the studies found in the literature), all

τ_{1 i}

values would be correlated with the standard deviation in measurements of anchor 1 (

σ_{1}^{2}

). Therefore, the resulting FIM for the measurement set

τ_{1} = \{τ_{12}, τ_{13}, τ_{14}\}

is as in Equation (A18).

F_{τ} = [\begin{matrix} s_{1} + s_{2} & s_{1} & s_{1} \\ s_{1} & s_{1} + s_{3} & s_{1} \\ s_{1} & s_{1} & s_{1} + s_{4} \end{matrix}]

(A18)

where

s_{i} = σ_{i}^{2}

Appendix C. Initial Filter Design

As discussed in Section 3, the IPS in our experimental setting is already equipped with an Extended Kalman Filter (EKF). The EKF developers note that the position estimates are affected by a non-uniform bias.

Our original idea was to use multiple filtering layers to enhance both precision and accuracy (see Figure A3). While the EKF is the first filter by default, the remaining ones may be applied in any order. The reason the debiasing filter (refer to Section 3.2) is applied last is that the bias values can change dramatically throughout the flying area, while the precision values do not (and have smaller magnitudes). Hence, we deemed it better to work on a more precise estimate to avoid selecting the wrong value of the bias. The proposed filtering process is defined as follows:

1.: Extended Kalman Filter
2.: Saturation (and artificial smoothing)
3.: Correction of position via fourth-order Adams–Moulton (AM4) method
4.: Debiasing filter

Figure A3. Proposed filtering process, consisting of four steps and returning estimates identified with their respective filter symbols:

\bar{x}

stands for saturated,

\tilde{x}

for dynamically corrected, and

\hat{x}

for debiased. Superscripts

^{b}

and

^{g}

refer to the body frame and inertial frame, respectively.

Figure A3. Proposed filtering process, consisting of four steps and returning estimates identified with their respective filter symbols:

\bar{x}

stands for saturated,

\tilde{x}

for dynamically corrected, and

\hat{x}

for debiased. Superscripts

^{b}

and

^{g}

refer to the body frame and inertial frame, respectively.

Unfortunately, precision was not noticeably enhanced. Therefore, only the EKF and the Debiasing Filter discussed in Section 3 were implemented for the experiments carried out in this paper (see Section 4 and Section 5). Nonetheless, the mathematical formulations of the other filtering layers are included below to support future work.

Appendix C.1. Saturation and Smoothing Filter

The saturation filter limits the eventual estimation overshoots of both velocity and position. Two types of velocity overshoots are considered, namely, a maximal overestimation and a maximum allowed time derivative. For example, each component of the estimated velocity vector is limited by a maximum velocity value that can be different in the three directions. For instance, the horizontal velocity components (in x and y direction) can be limited by the maximum cruise speed, while the vertical velocity (in z direction) can be at maximum two times the free-fall speed. Moreover, the time derivative of the velocity is limited by the expected maximum acceleration. For simplicity of analysis, the derivatives in all directions have been limited by the same amount, though this may not be the case in practice. The effect of the presented saturation filter on a casual sequence of velocity estimations is depicted in Figure A4. The filter acts analogously on a sequence of position estimates, i.e., Equations (A23)–(A25).

Figure A4. Representation of the presented saturation filter on velocity estimates in ith direction; blue shows the original EKF estimate and red shows the correction. Note that the measurements are discrete and represented by the peaks; the linear interpolation between measurements is only for visualisation purposes.

We start by defining the scalar variation of velocity (

^{g} Δ v^{(k)}

) as the modulus of the difference between the EKF-estimated velocity at the current time

^{(k)}

and the filtered velocity at the previous time step

^{(k - 1)}

^{g} Δ v^{(k)} = ∥Δ^{g} v^{(k)}∥ = ∥^{g} v^{(k)} -^{g} {\bar{v}}^{(k - 1)}∥

(A19)

This value should not exceed the maximum allowed speed variation (limited by the maximum acceleration,

a_{\max}

). Therefore, the inequality (A20) can be used to formulate the ceiling of the velocity variation (

Δ^{g} {\bar{v}}^{(k)}

) in Equation (A21):

^{g} Δ v^{(k)} ⩽ Δ t \cdot a_{\max}

(A20)

Δ^{g} {\bar{v}}^{(k)} = \min (Δ t \cdot a_{\max} \cdot {^{g} Δ v^{(k)}}^{- 1}, 1) \cdot Δ^{g} v^{(k)}

(A21)

At this point, every velocity component has to be limited by the maximum allowed speeds (

v_{i \max}

). The components of the final filtered velocity vectors (

^{g} {\bar{v}}^{(k)}

) are therefore defined as follows:

\begin{matrix} ^{g} {\bar{v}}_{i}^{(k)} = \min (^{g} {\bar{v}}_{i}^{(k - 1)} + Δ^{g} {\bar{v}}_{i}^{(k)}, v_{i \max}) \\ v_{\max} = [\begin{matrix} v_{1 \max} & v_{2 \max} & v_{3 \max} \end{matrix}] \end{matrix}

(A22)

The same filtering procedure can be applied to the position estimates, although in this case only the time derivative is limited. As before, the inequality (A20) can be used to define the variation of each

i^{th}

component (

i = 1, 2, 3

) of the filtered position

Δ {\bar{x}}^{(k)}

expressed in Equation (A24).

\begin{matrix} Δ x_{i}^{(k)} ⩽ Δ t \cdot v_{i \max} \end{matrix}

(A23)

where

Δ x^{(k)} = ∥Δ x^{(k)}∥ = ∥x^{(k)} - {\bar{x}}^{(k - 1)}∥

Δ {\bar{x}}_{i}^{(k)} = \min (Δ t \cdot v_{i \max} \cdot {Δ x^{(k)}}^{- 1}, 1) \cdot Δ x_{i}^{(k)}

(A24)

Finally the filtered estimation of the position can be defined as follows:

{\bar{x}}^{(k)} = {\bar{x}}^{(k - 1)} + Δ {\bar{x}}^{(k)}

(A25)

Embodied in this filtering step is the smoothing filter that aims to artificially reduce the oscillations of the measurements by averaging the history of

n + 1

measurements. The average is weighted such that the contribution of the last measurement (

{\bar{x}}^{(k)}

) is more important:

\begin{matrix} {\bar{x}}_{s}^{(k)} = {(1 - q)}^{n} \cdot {\bar{x}}^{(k - n)} + \sum_{n - 1}^{l = 0} q \cdot {(1 - q)}^{l} \cdot {\bar{x}}^{(k - l)} \\ q = q_{%} / 100 \end{matrix}

(A26)

where

q_{%}

is the percentage influence of the last measurement.

Appendix C.2. Adams–Moulton Filter

The main idea of Adams–Moulton dynamic prediction is to correct the EKF estimations, which consider only the measurements at the current time, using a prediction that takes into account the previous history of estimations. Hence, it is a multistep measurement filter. Usual multi-step predictor–corrector schemes consist of the combination of an explicit predictor, e.g., Adams–Bashforth (AB) of order n, with an implicit corrector, e.g., Adams–Moulton (AM) of order m, in order to take advantage of the prediction of the derivative function provided by AB(n). For the exception of a presented correction filter, AM can be used directly. In (A27), an AM scheme of fourth order is defined for predicting the position

{\bar{x}}_{AM 4}^{(k)}

via the use of filtered velocity estimates at the current and previous four time steps. A higher order AM scheme could be selected to consider a longer history of estimates.

\begin{matrix} {\bar{x}}_{AM 4}^{(k)} = {\bar{x}}^{(k - 1)} + \frac{Δ t}{720} \cdot (251 \cdot^{g} {\bar{v}}^{(k)} + 646 \cdot^{g} {\bar{v}}^{(k - 1)} \\ - 264 \cdot^{g} {\bar{v}}^{(k - 2)} + 106 \cdot^{g} {\bar{v}}^{(k - 3)} - 19 \cdot^{g} {\bar{v}}^{(k - 4)}) \end{matrix}

(A27)

One aspect that is considered by the EKF estimates and neglected by AM4 prediction is the effect that the control command input has on the state. For instance, the EKF takes into account the given thrust input. Intuitively, at low cruise speeds, the drone’s dynamics are very affected by control input; thus, the EKF estimate is more important, while at high cruise speed the drone dynamics at short period are mostly a direct effect of the stored momentum of the flying body, and move by inertia. Therefore, a vectorial weighting funcion

α

is defined here to ensure that the filtered correction of the estimated position follows this concept. The formulation of the newly filtered position

{\tilde{\bar{x}}}^{(k)}

is shown in (A28), while the shape of the ith weight component function (A29) is shown in Figure A5.

{\tilde{\bar{x}}}^{(k)} = α ⊙ {\bar{x}}^{(k)} + (1 - α) ⊙ {\bar{x}}_{AM 4}^{(k)}

(A28)

where ⊙ is the component-wise multiplication operation; thus, given two vectors

a

and

b

, the resulting vector components are

c_{i} = a_{i} \cdot b_{i}

, with

i = 1, 2, 3

α_{i} = exp (- \frac{1}{2} {(\frac{{\bar{a}}_{i}^{(k)}}{a_{i f}})}^{2}) (1 - α_{\min}) + α_{\min}

(A29)

Figure A5. Weight function for averaging EKF state estimations and AM4 predictions;

v_{i f}

is the flipping velocity,

v_{i \max}

is the maximum expected speed, and

α_{\min}

is a calibration parameter.

Figure A5. Weight function for averaging EKF state estimations and AM4 predictions;

v_{i f}

is the flipping velocity,

v_{i \max}

is the maximum expected speed, and

α_{\min}

is a calibration parameter.

References

Elsanhoury, M.; Makela, P.; Koljonen, J.; Valisuo, P.; Shamsuzzoha, A.; Mantere, T.; Elmusrati, M.; Kuusniemi, H. Precision Positioning for Smart Logistics Using Ultra-Wideband Technology-Based Indoor Navigation: A Review. IEEE Access 2022, 10, 44413–44445. [Google Scholar] [CrossRef]
Roy, P.; Chowdhury, C. A Survey of Machine Learning Techniques for Indoor Localization and Navigation Systems. J. Intell. Robot. Syst. 2021, 101, 63. [Google Scholar] [CrossRef]
Chen, D.; Neusypin, K.; Selezneva, M.; Mu, Z. New Algorithms for Autonomous Inertial Navigation Systems Correction with Precession Angle Sensors in Aircrafts. Sensors 2019, 19, 5016. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Rong, H.; Gao, Y.; Guan, L.; Zhang, Q.; Zhang, F.; Li, N. GAM-Based Mooring Alignment for SINS Based on An Improved CEEMD Denoising Method. Sensors 2019, 19, 3564. [Google Scholar] [CrossRef] [Green Version]
Widodo, S.; Shiigi, T.; Hayashi, N.; Kikuchi, H.; Yanagida, K.; Nakatsuchi, Y.; Ogawa, Y.; Kondo, N. Moving Object Localization Using Sound-Based Positioning System with Doppler Shift Compensation. Robotics 2013, 2, 36–53. [Google Scholar] [CrossRef] [Green Version]
Schott, D.J.; Saphala, A.; Fischer, G.; Xiong, W.; Gabbrielli, A.; Bordoy, J.; Höflinger, F.; Fischer, K.; Schindelhauer, C.; Rupitsch, S.J. Comparison of Direct Intersection and Sonogram Methods for Acoustic Indoor Localization of Persons. Sensors 2021, 21, 4465. [Google Scholar] [CrossRef] [PubMed]
Arbula, D.; Ljubic, S. Indoor Localization Based on Infrared Angle of Arrival Sensor Network. Sensors 2020, 20, 6278. [Google Scholar] [CrossRef]
Mahmoud, A.A.; Ahmad, Z.U.; Haas, O.C.; Rajbhandari, S. Precision indoor three-dimensional visible light positioning using receiver diversity and multi-layer perceptron neural network. IET Optoelectron. 2020, 14, 440–446. [Google Scholar] [CrossRef]
Alarifi, A.; Al-Salman, A.; Alsaleh, M.; Alnafessah, A.; Al-Hadhrami, S.; Al-Ammar, M.; Al-Khalifa, H. Ultra Wideband Indoor Positioning Technologies: Analysis and Recent Advances. Sensors 2016, 16, 707. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Tsang, Y.P.; Wu, C.H.; Ip, W.; Ho, G.; Tse, M. A Bluetooth-based Indoor Positioning System: A Simple and Rapid Approach. Annu. J. IIE 2015, 35, 11–26. [Google Scholar]
Zhao, X.; Xiao, Z.; Markham, A.; Trigoni, N.; Ren, Y. Does BTLE measure up against WiFi? A comparison of indoor location performance. In Proceedings of the European Wireless 2014; 20th European Wireless Conference, Barcelona, Spain, 14–16 May 2014; pp. 1–6. [Google Scholar]
Ezhumalai, B.; Song, M.; Park, K. An Efficient Indoor Positioning Method Based on Wi-Fi RSS Fingerprint and Classification Algorithm. Sensors 2021, 21, 3418. [Google Scholar] [CrossRef] [PubMed]
Abusara, A.; Hassan, M.S.; Ismail, M.H. Reduced-complexity fingerprinting in WLAN-based indoor positioning. Telecommun. Syst. 2016, 65, 407–417. [Google Scholar] [CrossRef]
Sahota, H.; Kumar, R. Sensor Localization Using Time of Arrival Measurements in a Multi-Media and Multi-Path Application of In-Situ Wireless Soil Sensing. Inventions 2021, 6, 16. [Google Scholar] [CrossRef]
Sakpere, W.; Oshin, M.A.; Mlitwa, N.B. A State-of-the-Art Survey of Indoor Positioning and Navigation Systems and Technologies. S. Afr. Comput. J. 2017, 29, 145–197. [Google Scholar] [CrossRef] [Green Version]
Hernandez, L.A.M.; Arteaga, S.P.; Perez, G.S.; Orozco, A.L.S.; Villalba, L.J.G. Outdoor Location of Mobile Devices Using Trilateration Algorithms for Emergency Services. IEEE Access 2019, 7, 52052–52059. [Google Scholar] [CrossRef]
Mosleh, M.F.; Zaiter, M.J.; Hashim, A.H. Position Estimation Using Trilateration based on ToA/RSS and AoA Measurement. J. Phys. Conf. Ser. 2021, 1773, 012002. [Google Scholar] [CrossRef]
Neirynck, D.; Luk, E.; McLaughlin, M. An alternative double-sided two-way ranging method. In Proceedings of the 2016 13th Workshop on Positioning, Navigation and Communications (WPNC), Bremen, Germany, 19–20 October 2016. [Google Scholar] [CrossRef]
Jamil, F.; Iqbal, N.; Ahmad, S.; Kim, D.H. Toward Accurate Position Estimation Using Learning to Prediction Algorithm in Indoor Navigation. Sensors 2020, 20, 4410. [Google Scholar] [CrossRef] [PubMed]
Mahida, P.; Shahrestani, S.; Cheung, H. Deep Learning-Based Positioning of Visually Impaired People in Indoor Environments. Sensors 2020, 20, 6238. [Google Scholar] [CrossRef] [PubMed]
Alraih, S.; Alhammadi, A.; Shayea, I.; Al-Samman, A.M. Improving accuracy in indoor localization system using fingerprinting technique. In Proceedings of the 2017 International Conference on Information and Communication Technology Convergence (ICTC), Jeju, Republic of Korea, 18–20 October 2017. [Google Scholar] [CrossRef]
Alhammadi, A.; Alraih, S.; Hashim, F.; Rasid, M.F.A. Robust 3D Indoor Positioning System Based on Radio Map Using Bayesian Network. In Proceedings of the 2019 IEEE 5th World Forum on Internet of Things (WF-IoT), Limerick, Ireland, 15–18 April 2019. [Google Scholar] [CrossRef]
Alhammadi, A.; Hashim, F.; Rasid, M.F.A.; Alraih, S. A three-dimensional pattern recognition localization system based on a Bayesian graphical model. Int. J. Distrib. Sens. Netw. 2020, 16, 155014771988489. [Google Scholar] [CrossRef]
European Union Agency for the Space Programme. Galileo Initial Services. 2021. Available online: https://www.euspa.europa.eu/european-space/galileo/services/initial-services (accessed on 18 August 2021).
Yeh, S.C.; Hsu, W.H.; Su, M.Y.; Chen, C.H.; Liu, K.H. A study on outdoor positioning technology using GPS and WiFi networks. In Proceedings of the 2009 International Conference on Networking, Sensing and Control, Okayama, Japan, 26–29 March 2009. [Google Scholar] [CrossRef]
Ghavami, M.; Michael, L.; Kohno, R. Ultra Wideband Signals and Systems in Communication Engineering; John Wiley & Sons: Hoboken, NJ, USA, 2007. [Google Scholar]
Innocente, M.S.; Grasso, P. Self-organising swarms of firefighting drones: Harnessing the power of collective intelligence in decentralised multi-robot systems. J. Comput. Sci. 2019, 34, 80–101. [Google Scholar] [CrossRef]
Bezas, K.; Tsoumanis, G.; Angelis, C.T.; Oikonomou, K. Coverage Path Planning and Point-of-Interest Detection Using Autonomous Drone Swarms. Sensors 2022, 22, 7551. [Google Scholar] [CrossRef] [PubMed]
Cramér, H. Mathematical Methods of Statistics (PMS-9); Princeton University Press: Princeton, NJ, USA, 1999; Volume 9. [Google Scholar]
Chen, C.S.; Chiu, Y.J.; Lee, C.T.; Lin, J.M. Calculation of Weighted Geometric Dilution of Precision. J. Appl. Math. 2013, 2013, 953048. [Google Scholar] [CrossRef] [Green Version]
Sieskul, B.; Kaiser, T. Cramer-Rao Bound for TOA Estimations in UWB Positioning Systems. In Proceedings of the 2005 IEEE International Conference on Ultra-Wideband, Zurich, Switzerland, 5–8 September 2005. [Google Scholar] [CrossRef]
Amigo, A.G.; Closas, P.; Mallat, A.; Vandendorpe, L. Cramér-Rao Bound analysis of UWB based Localization Approaches. In Proceedings of the 2014 IEEE International Conference on Ultra-WideBand (ICUWB), Paris, France, 1–3 September 2014. [Google Scholar] [CrossRef]
Zhang, J.; Kennedy, R.A.; Abhayapala, T.D. Cramér-Rao Lower Bounds for the Synchronization of UWB Signals. EURASIP J. Adv. Signal Process. 2005, 2005, 293649. [Google Scholar] [CrossRef] [Green Version]
Alhakim, R.; Raoof, K.; Simeu, E.; Serrestou, Y. Cramer–Rao lower bounds and maximum likelihood timing synchronization for dirty template UWB communications. Signal Image Video Process. 2011, 7, 741–757. [Google Scholar] [CrossRef]
D’Amico, A.A.; Mengali, U.; Taponecco, L. Cramer-Rao Bound for Clock Drift in UWB Ranging Systems. IEEE Wirel. Commun. Lett. 2013, 2, 591–594. [Google Scholar] [CrossRef]
Mallat, A.; Louveaux, J.; Vandendorpe, L. UWB based positioning: Cramer Rao bound for Angle of Arrival and comparison with Time of Arrival. In Proceedings of the 2006 Symposium on Communications and Vehicular Technology, Liege, Belgium, 23 November 2006. [Google Scholar] [CrossRef]
Silva, B.; Pang, Z.; Akerberg, J.; Neander, J.; Hancke, G. Experimental study of UWB-based high precision localization for industrial applications. In Proceedings of the 2014 IEEE International Conference on Ultra-WideBand (ICUWB), Paris, France, 1–3 September 2014. [Google Scholar] [CrossRef]
Grasso, P.; Innocente, M.S. Theoretical study of signal and geometrical properties of Two-dimensional UWB-based Indoor Positioning Systems using TDoA. In Proceedings of the 2020 6th International Conference on Mechatronics and Robotics Engineering (ICMRE), Barcelona, Spain, 12–15 February 2020. [Google Scholar] [CrossRef]
Compagnoni, M.; Notari, R.; Antonacci, F.; Sarti, A. A comprehensive analysis of the geometry of TDOA maps in localization problems. Inverse Probl. 2014, 30, 035004. [Google Scholar] [CrossRef] [Green Version]
Dulman, S.O.; Baggio, A.; Havinga, P.J.; Langendoen, K.G. A geometrical perspective on localization. In Proceedings of the First ACM International Workshop on Mobile Entity Localization and Tracking in GPS-Less Environments—MELT ’08, San Francisco, CA, USA, 19 September 2008. [Google Scholar] [CrossRef]
Park, K.; Kang, J.; Arjmandi, Z.; Shahbazi, M.; Sohn, G. Multilateration under Flip Ambiguity for UAV Positioning Using Ultrawide-Band. ISPRS Ann. Photogramm. Remote. Sens. Spat. Inf. Sci. 2020, V-1-2020, 317–323. [Google Scholar] [CrossRef]
Li, Y.L.; Shao, W.; You, L.; Wang, B.Z. An Improved PSO Algorithm and Its Application to UWB Antenna Design. IEEE Antennas Wirel. Propag. Lett. 2013, 12, 1236–1239. [Google Scholar] [CrossRef]
Lim, K.S.; Nagalingam, M.; Tan, C.P. Design and construction of microstrip UWB antenna with time domain analysis. Prog. Electromagn. Res. M 2008, 3, 153–164. [Google Scholar] [CrossRef] [Green Version]
Chahat, N.; Zhadobov, M.; Sauleau, R.; Ito, K. A Compact UWB Antenna for On-Body Applications. IEEE Trans. Antennas Propag. 2011, 59, 1123–1131. [Google Scholar] [CrossRef]
Chen, Z.N. UWB antennas: Design and application. In Proceedings of the 2007 6th International Conference on Information, Communications & Signal Processing, Singapore, 10–13 December 2007. [Google Scholar] [CrossRef]
Compagnoni, M.; Pini, A.; Canclini, A.; Bestagini, P.; Antonacci, F.; Tubaro, S.; Sarti, A. A Geometrical–Statistical Approach to Outlier Removal for TDOA Measurements. IEEE Trans. Signal Process. 2017, 65, 3960–3975. [Google Scholar] [CrossRef]
Compagnoni, M.; Notari, R. TDOA-based Localization in Two Dimensions: The Bifurcation Curve. Fundam. Inform. 2014, 135, 199–210. [Google Scholar] [CrossRef] [Green Version]
Kaune, R.; Horst, J.; Koch, W. Accuracy Analysis for TDOA Localization in Sensor Networks. In Proceedings of the 14th International Conference on Information Fusion, Chicago, IL, USA, 5–8 July 2011. [Google Scholar]
Kaune, R. Accuracy studies for TDOA and TOA localization. In Proceedings of the 2012 15th International Conference on Information Fusion, Singapore, 9–12 July 2012; pp. 408–415. [Google Scholar]
Grasso, P.; Innocente, M.S. Debiasing of position estimations of UWB-based TDoA indoor positioning system. In Proceedings of the UKRAS20 Conference: “Robots into the Real World” Proceedings, Lincoln, UK, 14–17 April 2020. [Google Scholar] [CrossRef]
Bitcraze. Loco Positioning System: TDOA Principles. 2020. Available online: http://www.bitcraze.io/documentation/repository/lps-node-firmware/2020.09/functional-areas/tdoa_principles/ (accessed on 18 August 2021).
Kay, S.M. Fundamentals of Statistical Processing; Prentice Hall: Hoboken, NJ, USA, 1993; Volume I. [Google Scholar]
Bitcraze. Loco Positioning System: TDOA2 vs. TDOA3. 2020. Available online: http://www.bitcraze.io/documentation/repository/lps-node-firmware/2020.09/functional-areas/tdoa2-vs-tdoa3/ (accessed on 18 August 2021).
DW1000 IEEE802.15.4-2011 UWB Transceiver—Datasheet v2.09. 2015.
DWM1000 IEEE 802.15.4-2011 UWB Transceiver Module—Datasheet v1.3. 2015.
Bitcraze. Getting Started with the Loco Positioning System. 2022. Available online: https://www.bitcraze.io/documentation/tutorials/getting-started-with-loco-positioning-system/ (accessed on 14 November 2022).
Mueller, M.W.; Hamer, M.; D’Andrea, R. Fusing ultra-wideband range measurements with accelerometers and rate gyroscopes for quadrocopter state estimation. In Proceedings of the 2015 IEEE International Conference on Robotics and Automation (ICRA), Seattle, WA, USA, 26–30 May 2015. [Google Scholar] [CrossRef]
Mueller, M.W.; Hehn, M.; D’Andrea, R. Covariance Correction Step for Kalman Filtering with an Attitude. J. Guid. Control. Dyn. 2017, 40, 2301–2306. [Google Scholar] [CrossRef]

Figure 1. Not-to-scale diagram of the

4 \times 4

m² 2D IPS being studied: (a) adjustable stands for the transmitting anchor antennas (A0–A3); (b) measurement points distributed every 50 cm in each direction; and (c) mobile stand for the object to be localised.

Figure 1. Not-to-scale diagram of the

4 \times 4

Figure 2. Dynamic experiment setup: mobile stand on rail; (a,b) beginning and end of rail; (h) optical obstacles (pins); (g) optical infrared sensor; (c) pcu, batteries, and motors; (d) Crazyflie 2.0; (e) direction of movement; (f) laser pointer. Optical sensor aligned with drone’s antenna.

Figure 3. (a) Original experimental radiation pattern sections on the

θ

ϕ_{1}

, and

ϕ_{2}

planes; (b) approximation procedure forcing identical values on intersections; and (c) radial projection of the approximated radiation pattern sections. These are used to reconstruct the 3D radiation pattern.

Figure 3. (a) Original experimental radiation pattern sections on the

θ

ϕ_{1}

, and

ϕ_{2}

Figure 4. Views of the reconstructed 3D radiation pattern of an anchor antenna.

Figure 5. Precision level sets and colourmaps for symmetric and random anchors. The magenta trapezoid is the convex hull of four anchors (modified from [38], with permission).

Figure 6. Bifurcation curves, bifurcation envelopes (green line), convex hull of anchors (magenta trapezoid, acceptable precision), and flyable area (yellow shade) for three and four anchors.

Figure 7. Expected position from the debiasing filter (

\hat{x}

) when applied to a measured posistion (

\bar{x}

) in 2D. The cloud of actual positions (

X_{i j}

) is constrained by the boundary

Ω

Figure 7. Expected position from the debiasing filter (

\hat{x}

) when applied to a measured posistion (

\bar{x}

) in 2D. The cloud of actual positions (

X_{i j}

) is constrained by the boundary

Ω

Figure 8. Representation of debiasing in 1D (left) and 2D (right).

Figure 9. Diagram of derivation of debiasing functions in x direction

^{x} β (x)

and y direction

^{y} β (x)

Figure 9. Diagram of derivation of debiasing functions in x direction

^{x} β (x)

and y direction

^{y} β (x)

Figure 10. Precision map of x component (

\pm^{x} σ

) and y component (

\pm^{y} σ

) of position.

Figure 10. Precision map of x component (

\pm^{x} σ

) and y component (

\pm^{y} σ

) of position.

Figure 11. Accuracy map of x component (

\pm^{x} b

) and y component (

\pm^{y} b

) of position.

Figure 11. Accuracy map of x component (

\pm^{x} b

) and y component (

\pm^{y} b

) of position.

Figure 12. Debiasing function for measured x component (

^{x} β

) and y component (

^{y} β

) of position.

Figure 12. Debiasing function for measured x component (

^{x} β

) and y component (

^{y} β

) of position.

Figure 13. Visualisation of rail and IPS position measurements, where

L_{r}

is the total length of the rail.

Figure 13. Visualisation of rail and IPS position measurements, where

L_{r}

is the total length of the rail.

Figure 14. Square path experiments (a) from [57] and (b) carried out in this paper. (a) Drone flying in auto-pilot along a desired square path (black) (from [57]). Both the EKF estimate (blue) and the actual position (red) are shifted from the desired path. (b) Proposed experiment following a fixed square path (black). The EKF + DF estimation (red) is expected to be more accurate than the EKF-only estimation (blue).

Figure 15. Absolute bias for x-direction measurements.

Figure 16. Absolute bias for y-direction measurements.

Figure 17. Dynamic experiment at average cruise velocity of 0.33 m/s with x spanning 0 m to 4 m at constant y = 2 m, where (

x, y

): position estimate with EKF-only; (

\hat{x}, \hat{y}

): debiased position; (

x_{ref}

): actual position on rail; and (

v_{x}

): estimated instantaneous velocity in x-direction.

Figure 17. Dynamic experiment at average cruise velocity of 0.33 m/s with x spanning 0 m to 4 m at constant y = 2 m, where (

x, y

): position estimate with EKF-only; (

\hat{x}, \hat{y}

): debiased position; (

x_{ref}

): actual position on rail; and (

v_{x}

): estimated instantaneous velocity in x-direction.

Figure 18. Square-path experiment results, with flying domain delimited by four anchors. The vehicle starts moving from

(0.5, 0.5)

following the positive x-axis direction. IPS-1 uses EKF only, while IPS-2 uses EKF + DF. Problematic regions of the path are highlighted in yellow. The overall experiment shows a clear improvement when incorporating the proposed DF.

Figure 18. Square-path experiment results, with flying domain delimited by four anchors. The vehicle starts moving from

(0.5, 0.5)

Table 1. Representative results of dynamic validation. The RMSEs of an IPS with and another without DF (IPS-2 and IPS-1, respectively) are compared to demonstrate the accuracy improvement of the former. The average performance difference is shown in columns

Δ x

and

Δ y

. The average cruise velocity is shown in the last column.

Δ x

and

Δ y

. The average cruise velocity is shown in the last column.

			RMSE_x,avg [cm]			RMSE_y,avg [cm]
dir.	$x$ [m]	y[m]	IPS-1	IPS-2	$Δ x$	IPS-1	IPS-2	$Δ y$	$v_{avg}$
hor.	[0, 4]	1	12.7	6.8	5.9	10.0	7.9	2.1	0.58
hor.	[0, 4]	2	12.0	8.1	3.9	6.7	4.3	2.4	0.44
hor.	[0, 4]	3	12.6	8.0	4.6	9.3	8.0	1.3	0.43
ver.	1	[4, 0]	15.6	10.3	5.4	9.4	6.8	2.7	0.58
ver.	2	[4, 0]	10.3	8.0	2.3	15.8	10.1	5.7	0.51
ver.	3	[4, 0]	11.4	9.4	2.0	15.3	12.2	3.1	0.42

Table 2. Square-path experiment results. The RMSEs of an IPS with and another without DF (IPS-2 and IPS-1, respectively) are compared to show the accuracy improvement. The ‘raw’ heading refers to full data stream and the ‘sel.’ heading refers to the undamaged data stream (i.e., no misbehaviour). Average improvement with DF is shown by

Δ

Δ

				RMSE_IPS-1 [cm]		RMSE_IPS-2 [cm]
Edge	dir.	x[m]	y[m]	Raw	sel.	Raw	sel.	$Δ$ [cm]
bot	hor.	[0.5, 3.5]	0.5	9.2	9.5	7.5	4.7	4.8
right	ver.	3.5	[0.5, 3.5]	12.6	12.5	9.0	8.3	4.2
top	hor.	[3.5, 0.5]	3.5	6.0	5.5	5.7	4.6	0.9
left	ver.	0.5	[3.5, 0.5]	15.8	15.2	8.8	6.7	8.5

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Grasso, P.; Innocente, M.S.; Tai, J.J.; Haas, O.; Dizqah, A.M. Analysis and Accuracy Improvement of UWB-TDoA-Based Indoor Positioning System. Sensors 2022, 22, 9136. https://doi.org/10.3390/s22239136

AMA Style

Grasso P, Innocente MS, Tai JJ, Haas O, Dizqah AM. Analysis and Accuracy Improvement of UWB-TDoA-Based Indoor Positioning System. Sensors. 2022; 22(23):9136. https://doi.org/10.3390/s22239136

Chicago/Turabian Style

Grasso, Paolo, Mauro S. Innocente, Jun Jet Tai, Olivier Haas, and Arash M. Dizqah. 2022. "Analysis and Accuracy Improvement of UWB-TDoA-Based Indoor Positioning System" Sensors 22, no. 23: 9136. https://doi.org/10.3390/s22239136

APA Style

Grasso, P., Innocente, M. S., Tai, J. J., Haas, O., & Dizqah, A. M. (2022). Analysis and Accuracy Improvement of UWB-TDoA-Based Indoor Positioning System. Sensors, 22(23), 9136. https://doi.org/10.3390/s22239136

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

Article Menu

Analysis and Accuracy Improvement of UWB-TDoA-Based Indoor Positioning System

Abstract

1. Introduction

2. Theoretical Study of Precision and Failure

2.1. IPS under Study

2.2. CRLB Analysis for Pseudo-Range Multilateration with Round-Robin Scheduling

2.2.1. Signal-to-Noise Ratio Formulation

2.2.2. Radiation Pattern of the DW1000 Anchor Antenna

2.2.3. Analytical Results of CRLB Analysis

2.3. Bifurcation Envelope Analysis

2.3.1. Bifurcation Curve

2.3.2. Bifurcation Envelope

3. Filter Design

3.1. Proposed Filter Design

3.2. Debiasing Filter

4. Design of Experiments

4.1. IPS Bias Map Generation

4.2. DF Calibration and Validation Setup

4.3. DF Validation under Dynamic Setup Conditions

4.4. Square Path Experiment Setup

5. Results and Discussion

5.1. Proof of Accuracy Improvement

5.2. Dynamic Validation of Debiasing

6. Conclusions and Future Work

Author Contributions

Funding

Conflicts of Interest

Abbreviations

Appendix A. Radial Basis Function Network Implementation

Appendix B. Reference CRLB Analysis

Appendix C. Initial Filter Design

Appendix C.1. Saturation and Smoothing Filter

Appendix C.2. Adams–Moulton Filter

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI