Efficient 3D trilateration algorithm for visible light positioning

This letter presents an efficient algorithm for estimating the three-dimensional (3D) location of a photodiode (PD) receiver via visible light positioning. It solely works on measured powers from different light-emitting diode (LED) sources and does not require any prior knowledge of the PD receiver height. It is found that four LEDs are required that are not on the same circle, in order to unambiguously determine the 3D location. The algorithm is optimized towards a minimized calculation time in view of real-time operation on energy-constrained lightweight and mobile devices such as drones.


Introduction
The advent of light-emitting diodes (LEDs) has sparked a large research interest in visible light positioning (VLP), whereby the location of a photodiode (PD) receiver is being estimated based on its received powers from different LEDs. For many applications, the PD is expected to maintain a fixed and known height, reducing the localization problem to a twodimensional (2D) problem. This is either solved via (manual or model-based) fingerprinting maps or via a classic trilateration method. For three-dimensional (3D) problems however, a fingerprinting approach with cm-level granularity becomes unfeasible due to the requirement of large memory and computation time to iterate over the 3D map. Alternatively, 3D trilateration can be applied, for which a plethora of algorithms is available. In VLP however, the distance cannot be estimated directly from the received power, as this estimation also depends on the height difference between the LED and the PD. Therefore, some workarounds have been sought, e.g. one can assume that the PD height is either known or estimated from on-board sensors, leading to a '2.5D solution'. In [1], a 3D algorithm is presented, using multiple receivers instead of multiple transmitters. In [2], a full 3D VLP algorithm is presented, estimating the 2D position based on the previous height estimate, followed by a height adjustment of this 2D location. However, this introduces errors as the best 2D estimate at the previous height will be different from the one at the new height. In [3], another full 3D approach is presented, but the method requires integral calculations and is limited to three sources, which will be shown further to not suffice for an unambiguous location estimation. The solution of [4] is also limited to three sources, and is mathematically relatively complex. Compared to available research, this letter presents a novel simple and error-free method for estimating the full 3D (x, y, z) location of a single PD, i.e. without any prior knowledge of the PD receiver height. It is shown that three LEDs, or even four LEDs on a square's corner, do not allow an unambiguous location estimation. Moreover, thanks to its low complexity and an additional fast search optimization, the method is suited for real-time 3D VLP operation.

VLP configuration
We assume a generic configuration with N LEDs, each horizontally placed at a fixed height h LED , denoted with LED i and with coordinates (x i , y i , h LED ), i=1 .. N. The Lambertian mode of LED i is denoted as m i and its transmitted power as P Ti . The PD, with area A, is also horizontally oriented and located at the unknown location (x, y, z) at an unknown distance d i from LED i . The angle of irradiance is defined as θ i and the angle of incidence as q ¢ , i as shown in figure 1(a). A setup with four LEDs in a 5×5×5 m 3 volume is considered; see figure 1(b). According to the channel model of [5], the power P Ri received from LED i , at (x, y, z) is given by: Inversely, given a measured power P Ri received at an unknown location, and knowing for horizontally oriented LEDs and PD (see figure 1(a)), the distance d i between the PD and LED i can be estimated as  d , i based on the power P Ri received from LED i (i=1 .. N): with Δh=h LED −z=the (unknown) height difference between LED i and the PD. As such,  d i cannot be directly determined from P Ri , due to the angle-dependent behavior of the transmitter: locations at different distances from the transmitter (and at different heights) might lead to the same P Ri measurement. Therefore, a classic trilateration algorithm cannot be straightforwardly applied. In the following, we will denote  d i as  d i (Δh), due to its height dependency.

Algorithm
The proposed algorithm is based on an iterative 2D trilateration combined with a nonlinear least squares (NLLS) minimization to find the most likely 3D location. Both parts will be explained hereafter.
2D trilateration-since the distances d i cannot be estimated directly, we can only use a 2D trilateration at an assumed Δh, for N LED sources. The squared horizontal distance between LED i and the PD is given by: After eliminating the quadratic terms in x 2 and y 2 by These linear equations in x and y can be written as For any given Δh, x y can then be estimated as using the Moore-Penrose pseudo-inverse of A: NLLS optimization-in this phase, the 2D trilateration process is repeated for each of the candidate heights h, varying between a minimum height h min and maximal height h max <h LED , with a height resolution R h . The estimated location is found at the minimum of C(h), which is defined as the average squared error between the estimated distancesd i (h) from (2) and the distances of the estimated location (x(h), y(h), h) from (5) [6]. It determines to what extent the assumed height h when performing the trilateration, was likely  - Algorithm 1. Summarizes the iterative trilateration approach with NLLS.

Algorithm performance
The algorithm performance is first evaluated for the LED configurations denoted as CFG1 to 4 in table 1. The PD height is set to 1 m (see figure 1(b)) with A=1 mm 2 . The LEDs are mounted at h LED =5 m, with a transmitted power of 10 W each, and m=1 (except LED 4 in CFG4). Candidate heights vary between 0 and 4 m, with a resolution R h =1 mm.
Without noise added to the P Ri values of (1), it is observed that with three LEDs (CFG1), there are always two possible 3D solutions for the receiver. Figure 2 shows the cost function C(h) for location (5.5, 7.5, 1), becoming zero at h=1 m, but also at h=2.5 m, leading to the location (6.25, 7.5, 2.5). Adding a fourth LED (with the same P T and m) is able to solve this ambiguity, provided that is not located on the circle formed by the first three LEDs (see figure 1(b)), which is due to the radiation pattern's geometrical properties. Figure 2 shows that C(h) for CFG1 is the same as for CFG2 (curves overlap). It can indeed be calculated with (1) that the received powers at the two locations are the same in CFG2, i.e. (P R1 , P R2 , P R3 , P R4 ) (nW)=(101, 28, 101, 28). As such, the traditional much-used square-shaped LED configuration with 4 LEDs on a circle (CFG2) is not preferred, as the fourth LED has no added value compared to CFG1 with 3 LEDs. Adding a LED that is not located on circles formed by other LED combinations (e.g. LED 4 in CFG3), eliminates the second minimum at h=2.5 m, as shown in figure 2. A second way of resolving the location ambiguity of CFG1, is mounting a fourth LED on the circle, but with a different mode m. Figure 2 indeed shows a unique minimum of C(h) at h=1 m in CFG4 (with LED 4 on the circle, but with m=3). However, an LED configuration with different m values for different LEDs is usually not preferred in view of lighting conditions. The aforementioned analysis is valid for any 3D PD location.
After showing the proper operation of the 3D trilateration algorithm, we now investigate the performance when Gaussian noise with a standard deviation of 1 nW is added to the idealized P Ri values of (1). Per height level, 1301 points inside the LED triangle formed by LED 1-2-3 are evaluated on a uniform grid of 10 cm (see figure 1(b)), in CFG1, CFG3, and CFG5 (using 3, 4, and 5 LEDs, respectively). For each of the locations, five simulations are done with random noise levels and using R h =1 cm. Figure 3 shows the cumulative distribution function (cdf) of the positioning errors for CFG1-3-5. At h=1 m and for CFG1 (3 LEDs), it can be seen that slightly less than 60% of the points is clustered around the correct location with errors below 50 cm, and the rest around the equally likely alternatives having errors up to 3 m. For CFG3 (4 LEDs) and CFG5 (5 LEDs), median errors decrease (from 18.7 cm to 7.3 and 4.8 cm respectively), a.o. because there is only one cost minimum (see figure 2). For h=3 m and CFG1, the error

Fast search optimization
Using MATLAB R2018a on a 2.1 GHz Intel i7 processor with 8 GB RAM, one location estimation for CFG3 takes 16.6 ms on average for a height resolution of 1 cm (401 candidate heights between 0 and 4 m). This indicates that the proposed method is able to estimate the 3D position in real-time, making it suitable for drone navigation using only a PD, for example. Since most 3D VLP applications will rely on constrained receiver modules that are mobile, wireless, and lightweight, it is crucial to minimize calculation time (and thus power consumption) of the algorithm as much as possible. Therefore, a fast search method is added to algorithm 1. Instead of iterating uniformly over all candidate heights between h min and h max (see line 2 of algorithm 1), we propose an adapted version of the onedimensional golden section search (GSS) algorithm [7]. This fast search method iteratively narrows the search interval inside which the extremum is located. Since GSS can only be applied to unimodal functions and C(h) can possibly be bimodal (see e.g. CFG3 in figure 2), we first coarsely determine the interval inside which the lowest minimum is located, by evaluating C(h) at candidate heights within the interval [h min , h max ], with a resolution of 30 cm, as it is observed that location ambiguity occurs for PD heights spaced more than 30 cm apart. capped at h min and h max , respectively. The GSS algorithm is terminated when the search interval is less than 1 cm wide, in order to get a same height resolution R h as in the original algorithm 1. The estimated PD height is set as the average of the final interval's bounds and (x, y) is estimated with the classic 2D trilateration (lines 4 and 5 of algorithm 1). Figure 3 shows that for a same configuration (i.e. CFG3 evaluated at h=1 m), this efficiency-optimized algorithm produces the same error cdf as the non-optimized version. However, besides maintaining the positioning performance, it reduces the average execution time over the receiver grid by 90%, from 16.6 to 1.7 ms per 3D position estimation. For real-life deployments, the execution time can possibly be further reduced, by narrowing the initial [h min , h max ] interval, i.e. with bounds determined by the previous PD height, the maximal vertical PD speed and the location update rate.

Conclusion
A novel method for full 3D VLP trilateration is presented, which does not require any prior knowledge of the receiver height. It combines 2D-trilateration with an NLLS approach and, for LEDs mounted in a plane, is shown to require 4 sources not on one circle, to unambiguously extract the 3D location. The algorithm's proper operation is shown for different configurations. Computation time is limited to approximately 17 ms and can be further reduced to less than 2 ms using a fast search algorithm, making it suitable for realtime operation. Future work includes an experimental evaluation of the algorithm, also in the presence of reflections, as well as optimal LED placement for a maximal accuracy.