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

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×4m2 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 ±3cm, and accuracy improved by about 15cm for static and 5cm for dynamic measurements, on average.


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 longrange 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.

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.

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 × 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 Sections 4 and 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.

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).
where h(x) is the range vector (e.g., distance between receiver and anchors), F τ is the covariance matrix of theτ measurements, and tr(·) 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).
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[·] stands for expectation.
Thus, the information matrix of the TDoA measurements set τ rr = {τ 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 σ 2 i .
Kaune et al. [48] suggest that the variance for a specific source is as in Equation (9): where SNR 0 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 0 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.

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 ): 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): (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 tdBm (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 tdBi (θ 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 TB w , where k B is the Boltzmann constant for radiation of a black body (≈1.38 × 10 −23 J/K).

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.  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.

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 ±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.

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.

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Ẽ(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Ẽ(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 }. (d) Four anchors randomly placed.

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 ijk , 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).

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.

Proposed Filter Design
In our experimental setting, the object to be positioned is a Crazyflie 2.0 nanoquadcopter 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.

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 asx, the measured posistion asx, and the cloud of actual positions as X ij , as shown in Equation (15) and Figure 7.
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.
Assuming R to be negligible, the real position can be obtained from the measurement 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 β(x) and y β(x).
x y x x X 12 For the purpose of the following derivation, continuous interpolating functions of the x and y biases, both along theî andĵ axes, have to be obtained (e.g., cubic splines). Therefore, from the biases of x measurements around any x ij position, two interpolating functions can be obtained, namely, x i b ij and x j b ij along theî andĵ axes, respectively. Likewise, ŷ i b ij and ŷ j b ij 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î direction for the bias of the measurement of the x-component of the position. Thus, the interpolated bias function x i b ij multiplied by the weighting probability distribution x γ ij 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 ij fall in the interval between (x − 3 x σ) ij and (x + 3 x σ) ij , where x σ ij is the standard deviation of the Gaussian distribution x γ ij . The same integral can be evaluated in theĵ 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). Figure 8. Representation of debiasing in 1D (left) and 2D (right).
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) 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 [x i,jȳ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.
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).
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.

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 ij ). 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 Figures 10 and 11.

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.

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): 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 (k) r ) in the inertial frame of reference can be obtained by projecting s (k) along the angle between the rail and the x-axis, θ r : if n (k) > n (k−1) then: Δx : x measured on rail : x measured by IPS : history of IPS pos. Figure 13. Visualisation of rail and IPS position measurements, where L r is the total length of the rail.
In addition, the velocity estimations with IPS-1 and IPS-2 are compared to the discrete average velocity on the rail: 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.

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.  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).

Proof of Accuracy Improvement
The results of the experiments in Section 4.2 are presented in this section along with the accompanying Figures 15 and 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 Figures 15 and 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.

Dynamic Validation of Debiasing
In this section, the results of the dynamic experiments in Sections 4.3 and 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 Tables 1 and 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) to (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).

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 ±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.

Abbreviations
The following abbreviations are used in this manuscript:

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.
The Gaussian activation function ( f ij ) is defined on every marker (X ij ) that was previously used for mapping purposes. Hence, every marker represents a node in the network.
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 The constant c 2,ij (θ) is a periodic function interpolating the values of c 2,mn on a neighbourhood of eight surrounding points (in 2D).
Therefore, after all the c 2,ij (θ) constants have been defined, the resulting RBFN (shown in Figures A1 and A2) is as provided in Equation (A7).

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: 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 τ ij 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]: 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τ at time step k can be expressed aŝ whereτ ij is the individual pseudo-range measurement using the two ToAs of the node to the ith and jth anchors, and τ ij is the real range difference. While the superscript (k) is omitted, it is implicitly inferred in further analyses. Please note thatτ is a column vector, not a matrix. Here,τ ij can have as many components as the binomial coefficient 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: where x is the index ofτ measurements and tr(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: 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: Hence, the total FIM for the location, considering that each standard deviation changes in space (σ i = ξ(x, y)), is as follows: where tr(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 τ is Using an efficient unbiased estimator, it is proven [52] that F τ is the measurement covariance matrix: 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 τ 1i values would be correlated with the standard deviation in measurements of anchor 1 (σ 2 1 ). Therefore, the resulting FIM for the measurement set τ 1 = {τ 12 , τ 13 , τ 14 } is as in Equation (A18).  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) : 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 (∆ gv(k) ) in Equation (A21): g ∆v (k) ∆t · a max (A20) ∆ gv(k) = min ∆t · a max · g ∆v (k) −1 , 1 · ∆ 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 ( gv(k) ) are therefore defined as follows: 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 ∆x (k) expressed in Equation (A24). ∆x Finally the filtered estimation of the position can be defined as follows: x (k) =x (k−1) + ∆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 (x (k) ) is more important: x (k) s = (1 − q) n ·x (k−n) + n−1 ∑ l=0 q · (1 − q) l ·x (k−l) q = q % /100 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 x (k) AM4 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.
x (k) AM4 =x (k−1) + ∆t 720 · 251 · gv(k) + 646 · gv(k−1) −264 · gv(k−2) + 106 · gv(k−3) − 19 · gv(k−4) (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 positionx (k) is shown in (A28), while the shape of the ith weight component function (A29) is shown in Figure A5.
where is the component-wise multiplication operation; thus, given two vectors a and b, the resulting vector components are c i = a i · b i , with i = 1, 2, 3.