Fusion of Time and Angle Measurements for Digital-Twin-Aided Probabilistic 3D Positioning

This paper considers the following digital twin-based technique to address the NLoS uplink positioning problem. Several BS perform AoA measurements on an uplink signal transmitted by a UE. Then, using a digital twin of the environment, ray tracing in these AoA directions is performed. The intersection of the rays is then the estimated position of the UE. This approach, illustrated in Figure 1, can be called “reverse/backward ray tracing”.

When dealing with noisy AoA measurements, we can follow this procedure: For each BS, we sample the AoD for the rays to be launched according to the statistical characteristics of the AoA measurement error. We refer to this method as the Monte Carlo approach. Thus, numerous rays are launched from each BS. The result of this ray-launching procedure is a collection of points in the xy plane at the elevation of the UE, indicating where the rays intersect the xy plane. The distribution fitting approach, explained in the following section, is then used to obtain the probability distributions for each BS. Finally, these distributions are to be multiplied to obtain the position probabilities.

The angle information can be obtained by processing the received signal via several signal-processing algorithms such as Delay-and-Sum [8], MUSIC [9], MVDR [11], ESPRIT [10]. Each algorithm may lead to a different error statistics.

A ray-tracing propagation model is considered to compute the path of the radio waves in the environment. It enables to compute the candidate trajectories between a transmitter and a receiver, even in NLoS situations.

In this paper, a ray has either an infinite length (i.e., case of infinite number of bounces) or a finite length. In the following sections, we propose to discretize the rays into points. Then, each point as an associated length corresponding to the length between the start of the ray and the point position. The term “length of a ray” refers to the last point of a discretized ray if the ray length is finite.

Finally, regarding the propagation of the rays, we consider only specular reflections. Please consult [2, Sec. II] for a discussion on this model assumption.

Let $X$ be a random variable representing the position of a UE. Let $\theta$ be a random variable representing the true AoA of the signal, and $Y$ a random variable representing the BS measurement(s) of the AoA. Regarding the notations, we use $p(\theta|y)$ for $p(\Theta=\theta|Y=y)$ and similarly $p(x|y)$ for $p(X=x|Y=y)$. The distribution $p(\theta|y)$ represents the statistics of the uplink AoA measurement error.

We let $n$ be the number of measurements such that $\mathbf{y}=[y_{1},y_{2},…,y_{i},…,y_{n}]$. The vector $\mathbf{y}$ comprises the measurement performed by all BS. For the sake of simplicity, we assume that there is one set of measurements per BS, and therefore $n$ BS. The measurement $y_{i}=[y_{i}^{1},y_{i}^{2},...]^{T}$ may contain measurements of several categories, such as the AoA and ToA (see below).

In the considered positioning problem, the goal is to approximate $p(x|\mathbf{y})$.

For the 2D positioning case, points are obtained as the intersection between the launched rays and the 2D plane of interest. The elevation of this plane is the (assumed to be) known elevation coordinate $z$. This is illustrated in Figure 2. There is one colour for the set of points obtained for each BS.

For the 3D positioning approach, we suggest discretizing the rays into points through sampling. This yields 3D points as shown in Figure 3. In this example, the angle of the rays are chosen in a Monte Carlo manner based on noisy AoA measurements and their error statistics.

For the 2D case, for each set of points obtained from each BS (one color in Figure 2) a parametric distribution is fitted. The distribution can for instance be a 2D Gaussian mixture model (GMM), as illustrated in Figure 4.

For the 3D case, the same approach can be considered. Instead of using a 2D probability distribution, a 3D distribution should be used. One can fit 3D GMM on the set of 3D points obtained for each BS. As in the 2D case, the position probabilities are then obtained by multiplying the obtained 3D probability density functions. The same fitting algorithm as in [2] (expectation-maximization) can be used to compute the parameters of the probability density functions. This is illustrated for the 3D case in Figure 5. Simulation results provided at the end of the paper indicate that the 3D GMM distribution is adapted for the considered set of 3D points.

In addition to AoA measurements, a BS may also have information related to a propagation time (PT) of the signal between the UE to localize and the given BS. This PT information can be obtained in several manners such as:

• •

  By computing a difference between a ToA and ToD of the signal to obtain the ToF. One important constraint when measuring the PT in this manner is to have a fine clock synchronization between the UE and the BS.

• •

  By considering the delay corresponding to the largest coefficients of the channel impulse response.

• •

  By analyzing the RSSI.

Let us denote a PT measurement at BS $i$ as $y^{1}_{i}$. This PT information can be mapped to a suited length of ray $L$ or equivalently to the length associated to a 3D point. One should simply use the formula: $L=y^{1}_{i}\cdot c$, where $c$ is the propagation speed of the signal, such as the speed of light.

Given a probability distribution on the PT, we can deduce a probability distribution $p(L|y^{1}_{i})$ on the length of the points or on the length of the rays (if it is finite). This distribution is analogous to the distribution $p(\theta|y^{2}_{i})$ (where $y^{2}_{i}$ is an AoA measurement).

In some cases, a relative propagation time (RPT) rather than a PT may be available. This RPT information can be obtained from a time TDoA of the signal between two BS. It is therefore computed as the difference of two measured ToA at two different BS. Similarly to the PT, the RPT can be mapped to the distance information by using the constant $c$.

The main advantage of the RPT is that no synchronization between the UE to localize and the BS is required. Only the BS should be synchronized. As the BS are usually higher-quality hardware and connected, this is a high advantage.

The RPT provides information about the length difference between two points corresponding to two BS (points of different colors in Figure 3), but not on the length of the points. An example of point-length difference is shown in Table I. In the example, there are 4 point lengths for the first BS (obtained from one or more rays launched from this BS) and 5 point lengths for the second BS (obtained from rays launched from this second BS). This results in $4\times 5=20$ possible pairs. Those values are to be compared with an expected point-length difference $\Delta L$ obtained from the RPT measurement of the two considered BS.

Given an error statistics on the RPT, one has access to the distribution $p(\Delta L|y_{i}^{1},y_{j}^{1})$, where $y_{i}^{1},y_{j}^{1}$ represents the ToA of the signal at BS 1 and BS 2, used to compute the TDoA and therefore the RPT and length difference. As a result, the points of Figure 3 can be discriminated based on to their length difference:

• •

  Pairs of points having a likely length difference should be kept.

• •

  Pair of points having unlikely length difference should be discarded.

We explain how this operation can be realized in the next section.

Each point associated to a ray has the same probability given an AoA measurement but not given a PT measurement. Hence, an infinite-length ray is not associated to a probability distribution on the PT $p(L|y)$, but rather the points of the subsampled rays. As a result, to sample the points according to $p(L,\theta|y_{i}^{1},y_{i}^{2})$, one could proceed as described in the following Algorithm 1 (Monte Carlo approach).

Alternatively, one can switch step 3 and 4: For a given ray, one can first sample a distance according to $p(L|y_{i}^{1})$ and keep only the point having this distance on the launched ray.

This process enables to obtain, for each BS $i$, a set of points sampled according to $p(L,\theta|y_{i}^{1},y_{i}^{2})$.

As explained in Subsection IV-B, one aims at selecting pairs of points having a likely length difference. One should therefore first collect the length difference associated to all pairs of points, sort them by length, and select pairs of points by Monte Carlo sampling. The following Algorithm 2 describes this process. It is similar to Algorithm 1, but where the RPT is used instead of the PT.

This process enables to obtain, for each BS $i$, a set of points sampled according to the probability distribution $p(\Delta L,\theta|y_{i}^{1},y_{j}^{1},y_{i}^{2})$.

Both approaches described in Subsection V-A and V-B enables to down-select the points of Figure 3, which are selected only based on AoA measurements. For instance, using PT measurements for the down-selection leads to the subset of points shown on Figure 7.

Subsequently, this subset of points is to be used for distribution fitting described in Section III-B. This yields the fitted distributions (one distribution per BS) illustrated in Figure 8, which are more relevant compared to ones in Figure 5.

As a reminder, the position-probability estimates of the UE are then obtained by multiplying the obtained fitted distributions.

For AoA, PT, and RPT measurements, a Gaussian error model is considered: $p(\theta|y^{2})\sim N(\theta,\sigma_{\eta}^{2})$, $p(L|y^{1})\sim N(L,\sigma_{\nu}^{2})$, and $p(\Delta L|y_{i}^{1},y_{j}^{1})\sim N(\Delta L,\sigma_{\nu}^{2})$.

We use the same simulation environment as in [1][2]. The scene is illustrated in Figure 6 and has the following dimensions: width 8m, length 18m, and height 2.5m. This environment is inspired from the recommendations of the 5G Alliance for Connected Industries and Automation (5G-ACIA) for indoor industrial scenario [12]. The 4 BS are located in the top corners of the scene (see Figure 6). The UE is randomly dropped in the considered environment (i.e, with a uniform distribution). Hence, there is randomness both in the UE position and the BS measurement errors.

The positioning error is denoted by $\epsilon=||\hat{x}-x||$, where $\hat{x}$ corresponds to the position with the highest probability according to distribution obtained as the multiplication of the fitted GMM models and $x$ is the ground-truth position.

Figure 3 and 7 show the points obtained from using only AoA measurements and both AoA and PT measurements, respectively. When both AoA and PT measurements are used, the obtained points concentrate into a reduced number of zones. With AoA measurements only, the points are distributed along the propagation paths. The set of possible positioning solutions is therefore reduced with the additional PT measurements. This phenomenon is shown in Figure 5 and 8.

In Figure 9, the cumulative distribution function (CDF) of the positioning error is presented for cases both without and with PT measurements. The set of considered standard deviation of AoA noise is $\{$0.25°, 0.5°, 0.75°, 1°$\}$ and the standard deviation of the PT noise is 50 cm. We observe that the fusion of AoA and PT measurements significantly improves the positioning accuracy. For example, with a standard deviation of the AoA noise of 1°, adding PT measurement can lower the positioning error of 250 cm at the 90th percentile.

Figure 9 also shows the effect of increasing the AoA noise. The positioning error increases as the AoA noise becomes more significant. For example, with $\sigma_{\eta}=$0.25°, the 90th percentile is below 28 cm. With 1° this number increases to 67 cm.

In Figure 10, we report the CDF of the positioning error when the AoA noise standard deviation is 1° and when the PT noise standard deviation ranges in $\{20,30,50,70,100\}$ cm. As expected, reducing the noise of the PT measurements improves positioning accuracy. However, the impact of the PT noise is significantly reduced compared to the AoA noise. This can be explained as follows: Using AoA measurements only, multiple candidate positions are determined as the intersections of the likely modes of the GMM distributions. The PT information then helps distinguish between likely and unlikely intersections. As the number of such intersections is limited and may be widely spaced, even with high variance in the PT information, the candidate intersections can still be effectively differentiated.

Finally, Figure 11 reports simulation results when RPT measurements, instead of PT measurements, are combined with AoA. The performance improvement is less than that of the PT with the same standard deviation. This is not unexpected as the RPT measurements offers less direct information than the PT measurements. However, it is important to note that, from a system perspective, the RPT information is easier to acquire than the PT information. Nevertheless, using these additional RPT measurements enables better positioning performance than with only AoA measurements when the AoA measurement noise is high (1°).