Abstract
Unlike many urban localization methods that return point-valued estimates, a set-valued representation enables robustness by ensuring that a continuum of possible positions obeys safety constraints. One strategy with the potential for set-valued estimation is GNSS-based shadow matching (SM) in which one uses a three-dimensional (3D) map to compute GNSS shadows (where line-of-sight is blocked). However, SM requires a point-valued grid for computational tractability, with accuracy limited by grid resolution. We propose zonotope shadow matching (ZSM) for set-valued 3D-map-aided GNSS localization. ZSM represents buildings and GNSS shadows using constrained zonotopes, a convex polytope representation that enables propagating set-valued estimates using fast vector concatenation operations. Starting from a coarse set-valued position, ZSM refines the estimate depending on the receiver being inside or outside each shadow as judged by received carrier-to-noise density. We demonstrate our algorithm’s performance using simulated experiments on a simple 3D example map and on a dense 3D map of San Francisco.
1 INTRODUCTION
GNSS-based localization is often unreliable in dense urban areas. As illustrated in Figure 1, direct line-of-sight (LOS) GNSS signals can be blocked or reflected by tall buildings (Hofmann-Wellenhof et al., 1992), creating non-line-of-sight (NLOS) and multipath effects, thereby lowering the number of visible GNSS satellites available for localization (Zhu et al., 2018). In particular, LOS satellite signals in the cross-street direction are more likely to be blocked or reflected by buildings than signals along the street. Consequently, positioning accuracy is degraded in urban areas, especially in the cross-street direction.
NLOS and multipath effects in GNSS satellites lead to high localization uncertainty in urban areas, especially along the cross-street direction. The pink ellipse depicts uncertainty in the receiver position estimate. This figure is adapted from Wang et al. (2013b)
Robustness to degraded accuracy is important for safety-critical GNSS applications such as autonomous driving and drone delivery. As detailed later in our discussion of related work, set-valued position estimates enable robustness guarantees by ensuring an entire set of possible positions lies outside of obstacles or inside user-specified bounds, for example. Unfortunately, generating such position estimates is often challenging, requiring approximation with a discrete collection of points. In this work, we propose a novel method of generating set-valued estimates without discretization by combining recent three-dimensional (3D)-map-aided techniques from the GNSS community with recent geometric set representations from the controls community.
1.1 Related Work
A variety of approaches exist to address the challenge of positioning accuracy by leveraging a 3D map of an urban environment and buildings. For example, in ray tracing, one can create discrete position candidates with a particle filter and find all the possible paths from each particle to the GNSS satellite, while considering a limit on the number of interactions with nearby buildings per each transmitted ray (Miura et al., 2013; Suzuki, 2016). Essentially, ray tracing estimates the reflection route of NLOS signals and uses this to correct the bias in GNSS pseudoranges. While there exist commercial implementations of this approach (Iland et al., 2018), ray tracing is computationally-expensive and, instead, requires being offloaded to the cloud. One strategy to reduce the computational complexity is to apply machine learning (ML)-based GNSS to evaluate a predefined grid of position candidates with a neural network by collectively processing different heatmaps of pseudorange residuals (van Diggelen & Wang, 2020). This enables real-time performance with server-based processing, but achieves onboard tractability by considering a trade-off between the cross-street positioning accuracy and the computations executed to estimate receiver position.
An alternative approach is shadow matching (SM; Groves, 2011; Groves et al., 2015; Wang et al., 2013b), which can avoid offloading computation to the cloud (Wang et al., 2013a, 2014). SM uses a 3D map to estimate GNSS shadows of buildings, which are regions of the map where direct LOS signals are blocked. A GNSS shadow is depicted in Figure 2. The received carrier-to-noise density (C/N0) of GNSS satellite signals enables us to determine if a receiver is inside or outside the GNSS shadow of each building. However, to be computationally tractable, SM requires precomputing building skylines for a predefined grid of position candidates. These skylines are compared against the azimuth and elevation of GNSS satellites at runtime (Groves et al., 2015;Wang et al., 2013a, 2014). Each point-valued position candidate is evaluated by comparing the received C/N0 with the precomputed building skyline to estimate a visibility score that can be used to identify the most likely position candidate (or candidates) and also a weighted empirical covariance for downstream use (Wang et al., 2013a).
Our proposed technique on set-valued shadow matching using zonotopes to leverage constrained zonotopes and efficiently represent buildings, shadows, and a set-valued receiver position estimate; note that the 2D group plane is shown as having volume only for illustration’s sake.
A common aspect of the above methods is that they represent position estimates in a discretized way. Since accuracy depends on the number of discrete points, these methods face a challenge in the case of robust urban localization, or ensuring that a position estimate obeys user-specified protection levels or safety bounds. That is, one can use a finer discretization to achieve a more accurate estimate at the expense of more computation.
To avoid the point-valued discretization trade-off, one can instead consider a set-valued approach to compute a continuum of state estimates given the initial set of states and known measurement bounds (Combettes & Civanlar, 1991; Scott et al., 2016; Shi et al., 2015; Shiryaev & Podivilova, 2015). This enables robustness by ensuring that the set lies within, e.g., user-specified safety bounds. Recently, there has been growing interest in set-valued representations of receiver states and measurements for robust localization and motion planning in urban environments (Bhamidipati & Gao, 2020; Kousik et al., 2019; Shetty & Gao, 2020). Unfortunately, these methods typically require assuming a set-valued representation of uncertain measurements by overapproximating a confidence level set of Gaussian distribution using a polytope (Bhamidipati & Gao, 2020; Shetty & Gao, 2020), for example. In other words, one must make an assumption about the underlying distribution of measurements. Therefore, it remains an open challenge to generate set-valued receiver position estimates for the distributions that are common in urban environments (e.g., multi-modal with disconnected components).
One strategy for modeling arbitrary set-valued state estimates is to use zonotopes. Zonotopes are convex, symmetrical polytopes that can propagate set-valued state estimates using fast vector concatenation operations. Zonotopes are constructed as Minkowski sums (Althoff & Dolan, 2014) of line segments in an arbitrary-dimensional Euclidean space. These objects are well known for representing reachable sets of dynamic systems (Althoff, 2015), which can enable formally-verified motion planning, fault detection, and navigation (Bhamidipati & Gao, 2020; Kousik et al., 2019; Shetty & Gao, 2020). These objects can also be extended to constrained zonotopes, which can represent arbitrary convex polytopes (Raghuraman & Koeln, 2022; Scott et al., 2016), avoiding the limitation of symmetry. While zonotopes are widely used in robotics for path planning and collision avoidance, they have not been applied in the field of GNSS localization for addressing multipath/NLOS effects.
1.2 Proposed Method
Our current work proposes a novel paradigm for 3D-map-aided GNSS localization based on GNSS shadow matching: set-based shadow matching termed as zonotope shadow matching (ZSM). This work is based on our recent ION GNSS+ 2021 conference paper (Bhamidipati et al., 2021). We formulate buildings and GNSS shadows (regions where LOS signals are blocked) using constrained zonotopes wherein we compute the GNSS shadow for each satellite/building pair via the intersection of building shadow volume with that of ground plane. An illustration of these set operations is shown in Figure 2. The ZSM algorithm begins with a coarse set-valued position estimate, which we iteratively refine using the GNSS shadow of each satellite/building pair in an urban map. We use C/N0 values to decide whether the receiver is inside or outside each GNSS shadow and intersect the shadow with our position estimate if inside the shadow, or else subtract the shadow from our position estimate. The final set-valued position estimate from our proposed ZSM algorithm can be incorporated into downstream polytopic estimation or verification methods without requiring assumptions about the distribution of the estimate.
A key feature of our proposed ZSM is that it avoids gridding and, instead, directly estimates a set-valued receiver position estimate by leveraging constrained zonotopes and two-dimensional (2D) polytopes (i.e., polygons). Furthermore, given that zonotope operations need fewer computations, our approach requires significantly lower computational power during the offline step of shadow matching, while during online processing, the required computational power is on-par to that of conventional SM and lesser storage than that of conventional SM is incurred. Another advantage is that our formulation estimates a non-Gaussian distribution of possible receiver positions (a continuum) in an urban environment, based on many position statistics such as worst-case error, multiple/distinct ambiguous modes, error bounds of each mode, can be evaluated in an exact manner.
Given the recent availability of high-fidelity, open-source 3D building maps (Miura et al., 2015; Wang et al., 2013b), we assume the inaccuracies in the building boundaries and road lanes and the effect of building materials on GNSS signals to be negligible. We also consider an ideal LOS/NLOS classifier (i.e., with true positive 100% [LOS detection] and true negative 100% [NLOS detection]). Our choice of an ideal LOS classifier is justified as it provides an uncorrupted platform to compare the performance of our proposed set-valued shadow matching (i.e., zonotope shadow matching) with conventional SM. Also note that, in this work, we do not explicitly account for multipath effects. While this work demonstrates great success in achieving the exact bounds of set-valued receiver positioning for various simulation scenarios (see Section 4), ZSM suffers from challenges similar to that of shadow matching, including brittleness to misclassification and inability to resolve multi-modal ambiguities without additional information. Note that addressing these other challenges, which are mentioned in a detailed manner in Groves et al. (2015), is beyond the scope of this manuscript. However, we have proposed extensions of this current work that account for brittleness to satellite misclassification in Neamati et al. (2022a) and for resolving multi-modal ambiguities in Neamati et al. (2022b). While the ZSM method proposed in this work cannot directly handle satellite misclassifications associated with the state-of-the-art, off-the-shelf classifiers (Haosheng et al., 2020; Haosheng et al., 2018), we can still implement this work in real-world settings by leveraging a subset of GNSS satellites, either LOS or NLOS, that exhibit high satellite classifier probabilities (see experiments in Section 4). Also, given that our proposed ZSM method is an anytime implemen-tation that works in a sequential manner, we can achieve real-time performance via strategic ordering of GNSS satellites (e.g., with an objective to maximize the dilution of precision) so the more useful satellites are processed first and the algo-rithm is terminated when required to output the current estimate of the set-valued receiver position.
1.3 Contributions and Paper Organization
Our key contributions are as follows:
We propose a novel set-valued approach to representing GNSS shadows using zonotopes. Offline, we represent buildings using constrained zonotopes. At runtime, we generate the GNSS shadow for each satellite/building pair by extending the building in the shadow direction from the satellite to the building, as shown in Figure 2. This operation is efficient due to our representation of buildings and GNSS shadows as constrained zonotopes.
We propose the ZSM algorithm. Given an initial set-valued receiver position estimate, we iteratively subtract or intersect GNSS shadows across satellites to refine the estimate depending on the associated LOS/NLOS characteristics as judged by received C/N0 values. To our knowledge, this is the first work utilizing set representations for 3D-map-aided GNSS.
We experimentally demonstrate that ZSM computes a set-valued receiver position estimate unlike conventional SM (Groves, 2011; Groves et al., 2015; Wang et al., 2013a, 2013b) in which point-valued position accuracy is limited by its precomputed grid resolution. We perform experiment simulations using (a) a simple 3D map comprised of two urban buildings and (b) a dense 3D building map of San Francisco. We validate ZSM performance in terms of computation load (for both offline and online computations), point-valued estimation error, and the position bounds (size of the final set-valued receiver position estimate) in along-street and cross-street directions.
The remainder of the paper is organized as follows. Section 2 introduces mathematical notations for the relevant set representations. Section 3 states our proposed ZSM algorithm. Section 4 presents various experimental results that validate our proposed ZSM in simulation. Finally, Section 5 provides concluding remarks.
2 PRELIMINARIES OF SET REPRESENTATIONS
We now introduce the mathematical notation and definitions used throughout this paper and present constrained zonotopes.
2.1 Mathematical Notation and Definitions
We denote the natural numbers as 
and n-dimensional Euclidean space as 
. Points and scalars are italic lowercase; sets and matrices/arrays are italic uppercase. Sets of sets are script uppercase. An n × n identity matrix is In. An n × n array of zeros is 0n, and a same-sized array of ones is 1n. Similarly, an n × m array of zeros is 0n×m, and a same-sized array of ones is 1n×m. An empty vector or array is [].
Let 
. The Minkowski (set) sum (Althoff & Dolan, 2014) is ⊕, defined as A ⊕ B = {a + b | a ∈ A, b ∈ B}. For example, the Minkowski sum of two (non-parallel) line segments in 
is a convex polygon, whereas the Minkowski sum of two (non-parallel) polygons in 
is a convex polyhedron. The convex hull is convhull(A) = {λa1 + (1 – λ)a2 | a1, a2 ∈ A}. A convex polytope 
can be constructed as the convex hull of a set of vertices 
, which we call its vertex representation.
2.2 Constrained Zonotopes
We now define constrained zonotopes (Scott et al., 2016), which we use to represent buildings and GNSS shadows as shown in Figure 2. Consider a center 
, a generator matrix 
, and linear constraints defined by 
and 
. Per Scott et al. (2016), we define a constrained zonotope 
as a set:
1
The columns of G are called generators. Also note that a zonotope is a constrained zonotope with no constraints defined, which we denote by 
. An example 3D zonotope (i.e., no constraints) is shown in Figure 3(a).
Subfigure (a) is a zonotope (grey volume) with no constraints and three generators; subfigure (b) shows a 2D zonotope (blue) and a constrained zonotope (green) with constraints shown in red on the left.
For another example illustration with 
, and b = 1, Figure 3(b) shows the 2D zonotope zono(c, G, [], []) in blue while the 2D constrained zonotope zono(c, G, A, b) is shown in green.
3 PROPOSED ZSM ALGORITHM
Our proposed ZSM algorithm is summarized as follows. First, we create a 3D shadow volume for each satellite/building pair by extending the building in the shadow direction from the satellite to the building. Note, the building, shadow direction, and shadow volume are all represented as constrained zonotopes. Second, we intersect the shadow volumes with the ground plane (on which we assume the receiver is positioned) to compute GNSS shadows as 2D constrained zonotopes. Third, by comparing received C/N0 values against a user-specified threshold, we either subtract (if LOS signal) or intersect (if NLOS signal) each GNSS shadow with our current set-valued receiver position estimate. Starting from an initial set of position estimates, the final position estimate set is obtained by iterating through all buildings and satellites.
We note that this is a snapshot algorithm, meaning that we make use of only the instantaneous data available at any given time instance. Furthermore, the order in which GNSS satellites and buildings are considered does not affect the final set-valued estimate, meaning that the method is parallelizable. However, we leave extensions to time-varying shadows, receiver motion models, and parallelization as future work. To proceed, we describe our model and assumptions of the urban environment and GNSS satellite signals, then detail how we compute shadows and, finally, present the ZSM algorithm.
3.1 System Model
3.1.1 Urban Environment Map
We assume access to a 3D building map of the urban environment in which we plan to perform positioning, as is done in related work (Adjrad & Groves, 2017; Iland et al., 2018; van Diggelen & Wang, 2020). We assume such a map consists of a ground plane and a collection of buildings.
We assume that the ground plane, denoted 
, is either a single 2D constrained zonotope (as in Figure 2), or a collection of 2D constrained zonotopes: 
, where each Gk is a 2D constrained zonotope and 
. Thus, our algorithm does not require the ground plane to be constant within an area of interest. We can account for the variations in the ground height by modeling the ground plane as a collection of triangles, wherein each triangle is converted to a constrained zonotope, and the ground plane is represented by a collection of constrained zonotopes. There exist many standard ways to extract ground height, of which one straightforward way is to extract the ground plane information from the 3D building map being utilized. The other is to leverage the open-source information on digital terrain models (DTMs), digital surface models (DSMs), or digital elevation models (DEMs; Chen et al., 2017). We also assume the ground plane is bounded (implied by ngrnd being finite). In particular, we only consider a specific, bounded urban area that could be potentially identified by a coarse position estimate to within a few kilometers of the localization error.
We denote the set of all buildings as 
, where 
is finite and 
denotes the i-th building. We further assume that each building can be represented as a union of constrained zonotopes:
2
where ni is the number of constrained zonotopes representing the building Bi. Note that any convex polytope is a constrained zonotope (Scott et al., 2016), so this is not a restrictive assumption. Also note that, since the union of constrained zonotopes cannot be represented as a constrained zonotope (because it is not necessarily a convex polytope), in practice, we represent each building as a list of constrained zonotopes (i.e., the union is treated implicitly).
3.1.2 Preprocessing Standard 3D Maps
Representing buildings as constrained zonotopes requires preprocessing a 3D map. Standard 3D maps generated by computer-aided design software are often represented by a triangulation, or a union of triangles. That is, for a building 
:
3
such that 
is finite and each 
is a vertex of the triangle. Assuming our map is in such a standard format, we preprocess a 3D map as follows to satisfy our assumption in Equation (2).
First, note that the convex hull of constrained zonotopes is given by Raghuraman and Koeln’s (2022) Theorem 5 as follows. Let 
and 
. Suppose that 
, and 
. Then:
4a
4b
4c
Now, notice that each vertex can be represented as a constrained zonotope, tl,j = zono(tl, j, [],[], []). So, we can represent Tl as a constrained zonotope by taking convhull (convhull ({tl,1 j, tl,2}), {tl,3}), meaning that we apply Equation (4) twice. Thus, we preprocess a 3D map by iterating over all of its triangles and representing each one as a constrained zonotope. In this case, ntrng = ni. Note that another series of convex hull operations can be performed to combine the constrained zonotopes (each corresponds to a triangle) to form a single constrained zonotope associated with each building, in which case ni = 1.
3.1.3 Receiver Assumptions
We assume that the receiver is on the 2D ground plane, which we find reduces the computational burden and is common in conventional SM (Groves et al., 2015). We plan to extend this method to 3D in future work.
To create a set-valued receiver position estimate, we consider an area of interest (AOI):
5
where Ak is a single ground plane-constrained zonotope. Since the ground plane is bounded, the AOI is bounded. While it is reasonable to treat the entire ground plane as the AOI, we make this distinction to enable our algorithm to take advantage of the non-shadow-based localization data obtained from various information sources. Note that the requirements for determining the size of the AOI remain the same as those for conventional SM (i.e., the AOI needs to cover the user’s true location). Some potential methods for this include the following, among which a few are explained in a more detailed manner in Groves et al. (2015): (a) use of position estimate and uncertainty bounds from the conventional GNSS ranging solution, which is based on weighted analysis of all GPS satellites or only the LOS ranging signals; (b) use of other sensors, such as camera-based geotagging, WiFi-based positioning, sensor fusion of GPS-vision or GPS-lidar (GPS light detection and ranging), multi-agent system with inter-ranging; (c) use of predicted position estimates and uncertainty from onboard state estimation techniques, such as Kalman filters, particle filters; and (d) use of a heuristic approach that involves either considering a very large AOI and later narrowing down the correct mode from among the multiple ambiguous modes of shadow matching identified or starting with an initial AOI and increasing the size by a scale factor until a viable set-estimate of receiver position is obtained. Furthermore, the user can choose to eliminate the segments in the AOI that lie within building footprints since the user is expected to be outdoors and not inside the building.
3.1.4 Satellites and Received C/N0
We denote the observed set of GNSS satellites by 
, where each 
is a satellite location and nsat is the total number. We obtain the received C/N0 values from the GNSS receiver; we denote these values by 
, one for each satellite. In this work, we assume an ideal LOS/NLOS classifier and, thus, do not need to account for the associated false alarm or missed-detection rates.
3.2 Computing Shadows
Now we detail how we compute shadows using operations on constrained zonotopes.
3.2.1 Overview and Assumptions
To create shadows, first, we consider the shadow direction associated with each satellite/building pair. This is a unit vector pointing in the direction from the satellite’s position to a point in the building. Then, we compute the 3D shadow volume for the satellite/building pair by extending the building in the shadow direction. Finally, we compute the 2D GNSS shadow for the satellite/building pair, which is a subset of the AOI in which the receiver does or does not lie depending on the received C / N0.
Note, we have made a key assumption that the shadow direction can be treated as a single direction, as opposed to a cone of possible directions, because the distance from the satellite to any building is much larger than the size of the building. However, for a more accurate model, or if the satellite is at a low elevation angle that would produce a large shadow, we can instead consider a cone of shadow directions created by the convex hull of sj and Bi using Equation (4). Since this convex hull strategy requires a more complicated presentation of the algorithm, one can discard low-elevation satellites while implementing ZSM.
3.2.2 Computing Shadow Directions
We compute the shadow direction ℓi,j as a unit vector from the satellite position sj to a point bi ∈ Bi. Recall from Equation (2) that each building is represented as a union of constrained zonotopes. Recall also that we assume the shadow direction is the same for any bi ∈ Bi. So, we choose bi as the mean of the vertices of all of the constrained zonotopes comprising Bi.
To get the vertices of Bi, we use the following procedure. Recall that Bi is a union of constrained zonotopes. Using techniques from Althoff (2015) and Matt (2021), one can convert a constrained zonotope 
to a vertex representation as follows. Suppose Z has m generators (i.e., 
). Notice from Equation (1) that Z is an affine map of the polytope 
and 
. That is, we can write Z = c + GP. To get a vertex representation of Z, we first apply the work of Matt (2021) to enumerate the vertices of P and get a set 
for which P = convhull(W). Then, the vertices of Z are vertices of the set 
, which we compute using Althoff (2015).
We summarize the above procedure with a function that applies to the union of constrained zonotopes comprising Bi. Let 
be a set of constrained zonotopes. Then:
6
is a vertex representation of the polytope (not necessarily convex) created by the union 
.
Now, consider the set 
, where Zl is as in Equation (2). We estimate the point bi as 
. Finally, we compute the shadow direction as:
7
Next, we use ℓi,j to compute a building’s shadow as a volume in 3D space.
3.2.3 Computing Shadow Volumes
To compute a shadow volume for satellite/building pair (i, j), we extend Bi in the shadow direction ℓi,j This extension is achieved by taking the Minkowski sum of the building and the shadow direction.
Since the building and shadow direction are constrained zonotopes, to proceed, we now define the Minkowski sum of constrained zonotopes. Consider the example zonotopes Z1 = zono(c1, G1, A1, b1) and 
. Suppose that 
, and 
. Note, k1 is the number of constraints for Z1 and, similarly, k2 for Z2, as per Equation (1). Per Scott et al. (2016), the Minkowski sum of constrained zonotopes is:
8
Note, this follows from the definition of constrained zonotopes.
Now we can compute the shadow volume. First, consider the zonotope:
9
where ϵ denotes a scaling factor that is greater than at least the height of the tallest building in the set 
(in practice, we use ϵ = 105 meters). In other words, Li,j is a zonotope representing a long line segment extending in the shadow direction ℓi,j Then, the shadow volume of building Bi with respect to satellite j is:
10
where Zl are the constrained zonotopes comprising Bi as in Equation (2). In practice, we represent the union of constrained zonotopes as a list of constrained zonotopes; so, we hold on to each:
11
where we have omitted the i,j indices from Vl for ease of notation.
3.2.4 Computing GNSS Shadows
Finally, to compute a 2D GNSS shadow for the satellite/building pair (i, j), we intersect the shadow volume with the 2D AOI.
Recall that the shadow volume and the AOI are represented by constrained zonotopes. Again let Z1 = zono (c1, G1, A1, b1) and 
with 
, and 
. Per Scott et al. (2016), the intersection of constrained zonotopes is a constrained zonotope:
12
Therefore, for the shadow volume Vl and area of interest 
, we can compute a shadow Si,j. as:
13a
13b
where we use Equation (10) to expand Vi,j into a union of constrained zonotopes. As mentioned earlier, we represent a union of constrained zonotopes as a list of constrained zonotopes. So, in practice, we consider each k-th shadow
14
which is a single constrained zonotope corresponding to each Ak. We omit the i, j, l indices for ease of notation.
3.3 Algorithm Details
We now describe how to initialize and perform ZSM, summarized in Algorithm 1. This is a modification of the conventional SM approach originally introduced by Groves (2011). In particular, we enable set-valued receiver position estimates. We emphasize again that this is a snapshot method, considering only the data available at a single time instance.
// zonotope shadow matching (Snapshot)
3.3.1 Algorithm Overview
We begin with the entire AOI 
as an initial 2D set-valued estimate of our receiver position (Line 2). If we have no localization information available besides the fact that we are in an urban area, then we set 
. Then, for each GNSS satellite (Line 3), we perform the following steps:
First, for each building (Line 5), we compute the 3D shadow volume as in Equation (10; Lines 6–9). Then, we intersect the shadow volume with the AOI to find a 2D GNSS shadow (Line 11), which we efficiently convert to a generic vertex representation (Line 12).
Second, we concatenate the vertex-represented regions of 2D GNSS shadows across all the buildings (Line 13).
Third, if the C / N0 value for the current satellite is below a user-specified threshold, then the current satellite is NLOS, so we intersect the concatenated 2D GNSS shadow with our current set-valued position estimate (Line 15 and Figure 6[c]); otherwise, the satellite is LOS, so we subtract the GNSS shadow from our set-valued position estimate (Line 17 and Figure 6[d]). We assume the initial AOI (Line 2) to contain the true receiver position and, hence, the intersection of GNSS shadow with our current set-valued position estimate is never an empty set.
3.3.2 Computational Considerations
We present Algorithm 1 without parallelization for ease of exposition. However, one can parallelize Algorithm 1 as follows. When iterating over each i-th building and each j-th satellite, one can create a separate position estimate Pi,j for each (i, f) pair, then intersect or subtract each Pi,j from the initial position estimate P in Line 2 by checking if the corresponding building/satellite pair is NLOS or LOS as in Lines 14–17. Note, we use the unparallelized version of Algorithm 1 in our numerical experiments.
Note that efficient numerical tools, such as the Multi-Parametric Toolbox (Herceg et al., 2013) and MATLAB’s polyshape, exist to perform polytope intersection and subtraction operations in 2D (Lines 11–14). However, to enable leveraging these 2D polygon tools, we require the constrained zonotope intersection operation as in Equation (12) to intersect the 3D building shadow zonotopes with the 2D ground plane zonotopes.
4 EXPERIMENTAL RESULTS
We experimentally validated our proposed ZSM method in simulation. Particularly, we chose a simulated platform to conduct exhaustive analysis of our algorithm’s performance. First, we performed a preliminary experiment to illustrate the utility of our constrained zonotope representation. Then, we evaluated ZSM in comparison to conventional SM (Groves, 2011; Wang et al., 2013a) for two simulation experiments.
Our ZSM implementation was in MATLAB. We used the open-source Continuous Reachability Analyzer (CORA) toolbox (Althoff, 2015) to represent constrained zonotopes. We used the polyshape tool to perform set operations on 2D GNSS shadows. We utilized received C / N0 from only GPS; however note that our algorithm is generalizable to other GNSS constellations as well. As explained earlier in Section 1, we consider an ideal LOS/NLOS classifier in this work and, thus, do not account for satellite misclassifications.
4.1 Preliminary Experiment: Minkowski Sum Evaluation
To illustrate the utility of using constrained zonotopes, we compared the speed of the Minkowski sum on 1,000 randomly-generated polytopes with up to 100 vertices, which we represented as both a standard vertex representation and as constrained zonotopes. Recall that a convex polytope 
can be represented as the convex hull of a set of vertices 
as in Section 3.1.2. We used the vertex representation because 3D urban maps are typically represented as a triangulation of vertices as noted in Section 3.1.2. In other words, this comparison considered an alternative to ZSM wherein we would apply Algorithm 1 directly to a 3D map without converting it to constrained zonotopes.
We performed the Minkowski sum of both representations with a line segment to mimic the process of computing a building shadow zonotope and compareed the average computation time as shown in Figure 4. We implemented the vertex representation with MPT (Herceg et al., 2013) and the constrained zonotope representation with CORA (Althoff, 2015). This experiment was performed on a laptop computer with an 8-core 2.4 GHz CPU and 32 GB RAM.
Constrained zonotopes improve Minkowski sum computation time by an order of magnitude over a standard vertex representation (each red plus denotes an outlier).
The constrained zonotope representation of this operation (executed once) was an order of magnitude faster on average (71.3 ms for MPT vs. 5.5 ms for CORA) and had a smaller standard deviation (3.3 ms for MPT vs. 1.1 ms with fewer outliers for CORA). Since this Minokwski sum was performed multiple times (specifically, 
times) in Algorithm 1, we see that the ZSM formulation enables a huge speed increase over using other polytope representations for SM. We also note that the Minkowski sum for a constrained zonotope with a line segment always resulted in one additional generator per Figure 4, whereas it is unclear how many additional vertices might be required for the vertex representation. So it is possible to preallocate memory for the constrained zonotope to further increase performance.
4.2 Comparison Method and Validation Metrics
We validated our proposed ZSM via two simulation experiments: one using a simple 3D map comprised of two buildings and the other using a dense 3D building map of San Francisco. We compared ZSM’s performance with that of conventional SM, which was implemented in Groves (2011) and Wang et al. (2013a, 2013b). All computations were performed on a laptop with a 2-core 2.5 GHz CPU and 8 GB RAM. In these simulation experiments, we modeled each building as a constrained zonotope (i.e., ni = 1 in Equation [2]).
4.2.1 Details of Conventional SM
We applied the conventional SM technique presented in Wang et al. (2013a) and explained earlier in Section 1. Offline, the method considers a prespecified uniform grid of position candidates and then performs a precomputation step by tracing evenly-distributed azimuth and elevation rays for each position candidate and each building in a 3D map to compute predicted satellite visibility. The number of discrete candidates in the uniform grid is determined by its grid size, which is defined as the distance between any two candidates either in cross-street or along-street directions. Note that the distance between any two position candidates in the along-street direction is the same as in cross-street, thus each grid represents a square.
At runtime (online computation), conventional SM computes a visibility score at each position candidate by comparing the received C / N0 (used for LOS/NLOS classification) for each satellite to that of predicted satellite visibility (from offline computation). Then, the position candidates with the highest visibility scores (can be more than one) are identified as the most likely receiver positions. The conventional method also computes the weighted empirical covariance (Martens et al., 2003) by analyzing the position candidates and their associated weights based on the normalized visibility scores.
To analyze ZSM performance for a given initial AOI, we precomputed the visibility map of conventional SM for various grid sizes. As mentioned earlier in Section 1, a finer discretization is attained with a lower grid size, thus increasing the position accuracy of the conventional SM technique, but at the expense of higher computation cost (both offline and online).
4.2.2 Validation Metrics
We applied four validation metrics:
Offline computation load: For ZSM, this is the time required for converting buildings from vertices to constrained zonotopes. For conventional SM, this is the time required for precomputing the visibility map at each position candidate, and, thus, the predicted satellite visibility.
Online computation load: For ZSM, this is the time required to run Algorithm 1. For conventional SM, this is the time required to compute the visibility scores at all position candidates and compute the most likely position candidates (based on highest visibility scores) and the weighted empirical covariance.
Point-valued estimation error in cross-street and along-street directions: For ZSM, this is the error in cross-street and along-street directions between the true receiver location and centroid of the final set-valued position estimate (multiple centroid values obtained when disjoint components are present). For conventional SM, this is the error in cross-street and along-street directions between the most likely position candidates (highest visibility score) and the true receiver location.
Bounds in cross-street and along-street directions: For ZSM, this is the width of the bounding box that encloses the set-valued position estimate. Note that we report the width of individual bounding boxes when disjoint components are present in the final set-valued position estimate. For conventional SM, this is twice the 3σ (three standard deviations) bound of the visibility-based weighted empirical covariance in the cross-street and along-street directions. In particular, we evaluated the weighted empirical covariance across all the position candidates based on their visibility scores.
4.3 Simulated Experiment #1: Simple 3D Building Map
We first describe a simulation experiment that considers a simple 3D map comprised of two urban buildings. We use the real-world GPS data collected from an Android phone and, then, emulate the received C / N0 values with simulated NLOS effects.
4.3.1 GPS Data Set and NLOS Emulation
We used the openly available Google Android data set (Fu et al., 2020) with the data collection setup shown in Figure 5(a). To allow control over LOS/NLOS signals for each satellite, we used the data collected in open-sky conditions via a Pixel 4 XL Modded smartphone, which logged the GPS data, and a high-fidelity GNSS-RTK/IMU setup, which served as a reference ground truth.
Subfigure (a) is the experiment test platform utilized to validate the proposed ZSM via publicly-available Google Android data sets in Fu et al. (2020). Subfigure (b) shows the skyplot for a particular receiver location with nine visible GPS satellites. The GPS satellites are represented by circles with LOS satellites in blue and NLOS (simulated effects) in dark yellow.
To design an ideal LOS/NLOS classifier, we emulated NLOS effects as follows. First, we overlayed the test surroundings with simulated buildings to mimic an urban setting. We then chose a time when nine satellites were LOS, as judged by C / N0 values above 38-dB-Hz, which is representative of open-sky conditions based on the empirical studies conducted in Hetet (2000) and Kuusniemi et al. (2004). We induced simulated NLOS effects in a GPS satellite when the LOS vector between its location and the true receiver position was obstructed by any simulated building. In particular, we attenuated the received C / N0 values for these identified NLOS satellites to be below the 38-dB-Hz threshold. Given that we are considering an ideal LOS/NLOS classifier that is based on binary thresholding against 38 dB-Hz, we do not account for the emulation and inclusion of multipath effects in this work.
4.3.2 Experiment Setup
We set up the experiment as follows, also shown in Figure 6(a). For the urban 3D map, we considered two buildings. The receiver 2D position was initialized at (0, 0) m in local map coordinates and the initial AOI was chosen as the entire length and width of a street between the two buildings, shown later in Figure 6(b). We observed a total of nine visible GPS satellites whose skyplot is shown in Figure 5(b). As mentioned earlier, we induced simulated errors in received C / N0 to represent LOS/NLOS characteristics with respect to the true receiver position and 3D map. In Figure 5(b), the NLOS satellites are indicated by dark yellow circle markers whereas LOS ones are blue. For conventional SM, we considered three grid sizes with the distances between the discrete candidates in cross-street and along-street directions as follows: 15 m, 10 m, and 5 m. The number of position candidates for the grid sizes of 15 m, 10 m, and 5 m, were 63, 124, and 427, respectively.
4.3.3 Results and Discussion
The first few steps (top-down view) of Algorithm 1 are shown in Figures 6(b)–6(f). In particular, we illustrate the satellites, buildings, initial AOI, GNSS shadows, and intersection/subtraction procedure. The full results are summarized in Table 1 and the final localization results for proposed ZSM and conventional SM are illustrated in Figure 7(a) and Figure 7(b). For each grid size of conventional SM reported in Table 1, we reported the top three most-likely position candidates (ranked based on the highest visibility score and the least point-valued estimation error with respect to ground truth).
Subfigure (a) is an illustration of the experiment setup, while subfigures (b)-(f) represent iterations of the proposed ZSM algorithm. The magenta area is the receiver position estimate and the blue areas are shadows.
Finally, after iterating through all buildings and satellites, the algorithm output is a polygon (2D polytope) on the ground plane representing the set-valued estimate of receiver positions (potentially with multiple disjoint components) given the current snapshot of GNSS signals (Line 18). An illustration of ZSM in simulation is shown later in Figure 6.
The ZSM position estimate in subfigure (a) has a comparable estimation error to the conventional SM results in subfigures (b)-(d). Note that the ZSM error does not depend on a predefined grid density.
Performance Comparison Results for ZSM vs Conventional SM on a Simulation Experiment With two Buildings Seen in Figure 6(a) and Nine GPS Satellites seen in Figure 5(b).
Figure 6(a) shows the simulated setup with two urban buildings and nine GPS satellites. The black circle represents the true receiver position, blue circles represent LOS satellites, and LOS signals are represented by a green dashed line. NLOS satellites are represented by dark yellow circles with NLOS signals represented by a red dashed line. Figure 6(b) shows a top-down view with the initial AOI in magenta. Figures 6(c)–6(f) show the GNSS shadow for each satellite or building pair in cyan with the resulting receiver position estimate after the set intersection or difference in magenta. While we only illustrate the procedure for PRNS 29, 25, 26, and 12, the final result from all nine satellites is shown in Figure 7. Note that the satellites are plotted near the receiver for visualization only; the shadow directions as in Equation (7) were computed using the satellites’ actual positions. Given the received C / N0 and 3D map, we demonstrated that the proposed ZSM computes a set-valued receiver position estimate.
Figure 7(a) shows the ZSM receiver position estimate in magenta. Figures 7(b)–7(d) show the conventional SM with grid sizes of 5 m, 10 m, and 15 m, respectively, where the position candidates are color-coded using a jet colormap with blue indicating a low visibility score of 2 and red representing a high value of 9 (number of available GPS satellites). For a lower online computation time of 0.39 s, ZSM demonstrates a comparable point-valued estimation error and a smaller width than the densest conventional SM case in both the cross-street and along-street directions.
The results from the three cases of the conventional SM technique were used for validating the improved performance of our proposed ZSM algorithm. We saw a comparable point-valued estimation accuracy of ZSM (i.e., 3.46 m and 16.05 m in cross-street and along-street directions, respectively) as that of conventional SM with grid sizes of 5 m, 10 m, and 15 m, whose accuracies are reported in Table 1. Furthermore, unlike all the conventional SM cases, we showed that ZSM achieved a set-valued position estimate with an accuracy in width of only 17.87 m in the cross-street direction and 50.11 m in the along-street direction. We also found that the offline and online computations of ZSM only incurred an average runtime of 1.7 × 104 s and 0.39 s, respectively. Importantly, without the need for gridding, the ZSM result produced this point-valued positioning accuracy while returning a set-valued estimate, which contained the true receiver position. We anticipate that the result might be more accurate if we leveraged more GPS satellites and buildings that would provide additional geometry constraints in the shadow matching, thus reducing the centroid error and bounds of the disjoint component (ambiguous mode) containing the true receiver location, if any. This forms the premise of our next simulation experiment.
4.4 Simulated Experiment #2: Using 3D Building Map of San Francisco
We performed our next simulation experiment using a publicly-available 3D building map of San Francisco. For this dense 3D urban map, we validated our algorithm’s performance compared to that of the conventional SM technique using emulated GPS data. We also performed a sensitivity analysis of proposed ZSM by varying the number of buildings being considered.
4.4.1 GPS Data Set and NLOS Emulation
We generated GPS data using the simulation pipeline shown in Figure 8. Considering a desired static position (true receiver location), ephemeris file, and start time as inputs, we utilized a C++ language-based software-defined GPS simulator known as GPS-SIM-SDR (Bhamidipati et al., 2019; Ebinuma, 2018) to simulate the raw GPS samples, which are indicative of open-sky conditions. Later, we performed acquisition and tracking using a MATLAB-based software-defined radio known as SoftGNSS (Rojas, 2011) to generate received C / N0 values. Similar to NLOS emulation described in Section 4.3.1, we induced NLOS effects in C / N0 values based on the simulated buildings in San Francisco and true receiver location.
Our GPS software simulation pipeline uses two publicly-available software: GPS-SDR-SIM (Ebinuma, 2018) and SoftGNSS (Rojas, 2011). For a given true receiver location, we simulated the C/N0 values and, thereafter, induced NLOS effects using a 3D building map of San Francisco.
4.4.2 Experiment Setup
We utilized a publicly available 3D building map of San Francisco (György, 2018). As explained in Equation (3) of Section 3.1.2, standard 3D building maps are represented as union of triangles, with each triangle comprised of three vertices. An illustration of the vertex representation for the 3D building map of San Francisco is shown in Figure 9(a), wherein we isolated the vertices of each building using the theory of connected graphs (König, 1990) and color coded them independently. We then converted the entire 3D map from vertex representation to constrained zonotope representation, which is visualized in Figure 9(b). Note that the prominent landmarks in San Francisco, such as the Transamerica pyramid and Salesforce tower, can be easily identified in Figure 9(b). The storage space required for building-constrained zonotopes is quite less. For instance, the .mat file size for the eight buildings discussed in Section 4.4.4 is 5 KB while, for 14 buildings, it is 8 KB, and for 20 buildings, it is 11 KB. Similarly, the entire San Francisco map (shown in Figure 9[b]) is 345 KB. While this is acceptable for many practical applications, based on the user platform specifications, if further reduction in file size is required, one could employ alternate techniques such as storing a reduced order representation of building-constrained zonotopes or storing a smaller section of the map at the original order of the building-constrained zonotopes. Note that the experiments conducted in this section were fully simulated and, therefore, inaccuracies of this 3D building map compared to that of the San Francisco city are not relevant. This is because we simulated the true receiver location and classified the satellites as LOS/NLOS relative to this open-source 3D building map.
Preprocessing 3D building map of San Francisco: (a) shows the vertex representation, wherein the vertices of each building are separately color-coded; and (b) represents each building as a constrained zonotope, which is given as input to our proposed ZSM algorithm.
The true receiver position was initialized at (0, −18) m in local map coordinates. The initial AOI was chosen to lie within a size of 120 m × 120 m while excluding the regions that were within the building footprints. In particular, we chose the size of the AOI in a heuristic manner by penalizing the position solution and uncertainty bounds (3σ) estimated by the MATLAB-based SoftGNSS (Rojas, 2011), which was explained earlier in Section 4.4.1. A top-down view of this illustration in shown in Figure 10(a) and a relevant section of the 3D San Francisco map comprised of eight buildings (of which one of them was the Salesforce tower) is shown in Figure 10(b). We simulated 14 GPS satellites whose skyplot is shown in Figure 10(c), with LOS satellites indicated by blue circle markers and NLOS by dark yellow. For comparison with the conventional SM technique (explained in Section 4.2), we considered four grid sizes of 30 m, 10 m, 5 m, and 3 m. The number of position candidates for the grid sizes of 30 m, 10 m, 5 m, and 3 m, were 16, 97, 346, and 902, respectively.
Our second simulation experiment using a 3D map of San Francisco: Subfigures (a) and (b) show the initial AOI and surrounding buildings, while (c) shows the skyplot for the true receiver position with LOS satellites in blue and simulated NLOS in dark yellow.
4.4.3 Comparison Results with Conventional SM
Table 2 reports comparison statistics between the proposed ZSM with different satellite subsets and conventional SM with varying grid sizes in terms of offline and online computation load, point-valued estimation error with respect to true receiver location, and the associated position bounds. A visualization of our ZSM with all available 14 satellites, a random subset of 10 satellites, a random subset of 8 satellites, and satellites with elevation > 20° are shown in Figures 11(a), 11(b), 11(c), and 11(d), respectively. From among those illustrated earlier in Figure 10(c), the satellites with PRNs 1, 7, 14, and 15 were eliminated for the 10 random satellite case, while the satellites with PRNs 5, 7, 8, 9, 17, and 28 were omitted for the 8 random satellite case. Also, for the case with elevation > 20°, only the GPS satellites with PRNs 7, 13, 14, 17, 28, and 30 were considered. Intuitively, the reduced number of GPS satellites mimic the cases when a non-ideal LOS/NLOS classifier is in a simulation setting, wherein the users can execute the proposed ZSM by leveraging only the GPS satellites exhibiting a high probability of being either LOS or NLOS while omitting the less certain ones. The detailed steps of convergence for our proposed ZSM algorithm for all available 14 satellites can be viewed here: https://youtu.be/anIh4hd3ikw. Similarly, the illustrations of conventional SM with grid sizes of 3m, 5 m, 10 m, and 30 m are shown in Figure 12(a), 12(b), 12(c), and 12(d), respectively.
Performance Comparison Between ZSM and Conventional SM Based on a 3D Map of San Francisco
Our proposed ZSM algorithm estimates the receiver position with high precision (i.e., the area of the receiver estimate set is small), but produces a bimodal output in this example. The mode containing the true receiver position is not majorly affected by the consideration of different subsets of the available GPS satellites.
The results of conventional SM strongly depend upon grid size (a denser grid is more accurate). However, in this case, the uncertainty bounds in the conventional SM are over 100 m, indicating that ZSM is much more precise.
We validated that the proposed ZSM using all available 14 satellites successfully estimated two disjoint sets that showed a resemblance with the most-likely position candidates estimated using conventional SM of different grid sizes. For the conventional SM technique with a grid size of 30 m, we observed one most-likely candidate (highest visibility score) that showcased a point-valued accuracy of 12 m and 60 m in the cross-street and along-street directions, respectively. By implementing the conventional SM technique for a grid size of 10 m, we identified two most-likely candidates, both of which exhibited point-valued estimation error in a cross-street direction of 2 m. Similarly, conventional SM with grid sizes of 5 m and 3 m identified 1 – 3 grid points with a high visibility score of 14 in the near-vicinity of each of the two disjoint sets estimated via proposed ZSM. While increasing the grid resolution reduces the estimated bounds of conventional SM, we observed that these bound values were more than 100 m across all grid sizes in both along-street and cross-street directions. These large uncertainty bounds from conventional SM can be majorly attributed to the use of Gaussian distribution for approximating the non-linear, multi-modal shadow matching distribution.
In contrast, across all satellite subsets, our proposed ZSM successfully estimated the set containing the true receiver location as one of the outputs, without requiring discretizing. Particularly, when all 14 available satellites were utilized, ZSM demonstrated a higher point-value accuracy of 1.2 m in cross-street and 8.1 m in along-street direction, while achieving a small bound of only 5.8 m in cross-street and 16.2 m in along-street directions. We observed that, as the number of GPS satellites utilized in our proposed ZSM increased, the centroid error and the bounds of the mode that contained the true receiver location marginally decreased. Similarly, as the number of utilized GPS satellites increases, the offline computation load remains the same, but the online computation load marginally increases.
An interesting thing to note here is that, while ZSM estimates multiple disjoint components (modes), the ones that do not contain the true receiver location can be ruled out by fusing proposed ZSM with other information sources, say inertial measurement units (IMUs) or GPS pseudoranges. In particular, these additional sources could provide temporal information about the vehicle motion that could be used to perform a consistency check (i.e., compute the predicted/expected modes from IMU/GPS pseudoranges), and compare it with the set-valued estimate from zonotope shadow matching to eliminate incorrect modes over time. For instance, IMU data can be utilized to perform consistency checks over multiple time instants in the relative domain, itself, while considering the IMU biases and drifts to be constant over a relatively shorter time window. Similarly, GPS pseudoranges from satellites identified as LOS with high certainty can provide useful information, particularly in resolving the along-street component, thereby narrowing down the correct mode. Also, note that parallelization can further reduce the computational load of both the proposed ZSM and conventional SM.
4.4.4 Sensitivity Analysis of ZSM
We performed a sensitivity analysis of ZSM to evaluate the increase in offline and online computation loads as the level-of-detail for each building increased or alternatively the number of buildings increased. As explained in Section 3.1.2, standard 3D maps can be represented as a union of triangles, wherein each triangle is converted into a constrained zonotope. Therefore, in theory, the computational load incurred by increasing the number of buildings (with sparser level of detail for each building) is on the same order of magnitude as increasing the level of detail for each building, given a lower number of buildings. This is justified by comparing Equations (2) and (3) as follows: (a) given ntrng triangles and ni = ntrng, the total number of constrained zonotopes representing the buildings are ntrng (each building is sparser in the level of detail); and (b) given another ntrng triangles with ni < ntrng, the total number of constrained zonotopes representing the buildings would be ni (each building encompasses more level of detail). Given that the open-source 3D building map (György, 2018) is not dense enough to conduct a sensitivity analysis on level of detail explicitly, we analyzed the performance of the proposed ZSM as the number of buildings are varied. We also quantified our ZSM’s performance by analyzing the following: (a) centroid error and bounds associated with the disjoint component (ambiguous mode) containing the true receiver location; and (b) number of ambiguous modes and their variation as the number of buildings increases.
As shown in Table 3, we analyzed our ZSM’s performance by increasing the number of buildings from eight (discussed in Section 4.4.3) to 14 and 20. For all three cases, we considered the same initial AOI, which is indicated by the black box shown in Figure 13. Figure 13 illustrates the final set-valued position estimate obtained for the case with 20 buildings. We validated that our proposed ZSM technique was able to successfully detect all the disjoint components (ambiguous modes) and their exact bounds, which can vary based on the number, density, and placement of buildings considered. From Table 3, we observe that the number of disjoint components in the final set-valued position estimate increases as the number of buildings increases, wherein the number of disjoint components is two when eight buildings are considered, three disjoint components for 14 buildings, and five disjoint components when 20 buildings are considered. As explained earlier, disjoint components not containing the true receiver location can be omitted by fusing proposed ZSM with other sources, such as IMUs and GPS pseudoranges. We demonstrated that our proposed ZSM successfully maintained a high estimation accuracy and a small position bound in both cross-street and along-street directions. Another interesting point to note is that the eight-building and 14-satellite configuration already provides sufficient geometric diversity and, thus, increasing the number of buildings to 14 and 20 does not further reduce the centroid error and bound of the ambiguous mode containing the true receiver location.
ZSM produces a set-valued estimate with many disjoint components (magenta) when considering many buildings (pink) and a large AOI. However, each disjoint component is small, so fusing a ZSM estimate with other estimates to eliminate multi-modality can produce an accurate receiver position estimate.
ZSM Estimation Accuracy and Computation Time vs. Building Density
Based on the results, we observed that the offline computation load (even with no parallelization as of now) demonstrated easy scalability with a number of buildings (i.e., no significant increase in computation time). While the online computation load increases with the number of buildings, note that there is an extensive scope for parallelization in the extraction of GPS shadows that can drastically reduce the computation time. For instance, one can parallelize the processing of analyzing each pair of constrained zonotope (associated with 3D buildings) and satellite to extract the GPS shadow, which is currently executed in a sequential manner as explained in Algorithm 1. However, note that parallelization incurs an increase in the computational resources required.
4.4.5 Performance Analysis of ZSM: Different Location
We evaluated our ZSM’s performance in a different section of the 3D building map of San Francisco (i.e., near an intersection, where the distinction between the cross-street and along-street directions is not obvious). Given this, we labeled the directions as X and Y in the local map coordinates as shown in Figure 14(a), which comprised of 17 urban buildings in total. The true receiver position was initialized at (−20, −30) m in local map coordinates and the initial AOI (black box) was chosen to lie within a size of 212.1 m × 212.1 m while excluding the regions that lay within the building footprints. We considered 14 GPS satellites, which are same as those considered in earlier sections and whose skyplot is shown in Figure 14(b) with LOS satellites (a total of four) indicated by blue circle markers and NLOS by dark yellow. A visualization of our ZSM with all available 14 satellites and satellites with elevation > 20° are shown in Figures 14(c) and 14(d), respectively. Among all the GPS satellites shown in Figure 14(b), for the case with elevation > 20°, only the relevant six satellites with PRNs 7, 13, 14, 17, 28, and 30 were considered.
Our third simulation experiment in a different location using the 3D map of San Francisco; (a) shows a relevant section of 3D map near the Salesforce tower wherein the distinction between cross-street and along-street is not well-defined; (b) shows the skyplot with respect to the true receiver location with 14 visible GPS satellite, the GPS satellites are represented by circles with LOS satellites in blue and NLOS [simulated effects] in dark yellow; (c) shows the final set-valued estimate of the receiver position estimated by our ZSM when all 14 satellites are considered; and (d) shows this set-valued estimate using satellites with elevation > 20°. We demonstrate that in both cases, our ZSM successfully detects the mode containing the true receiver location, while exhibiting a higher centroid error and bounds for the reduced six satellite case as compared to the full 14 satellite case.
We demonstrated that our proposed ZSM estimated all the disjoint components (ambiguous modes). For the case in which all the GPS satellites are considered, three disjoint ones were identified with an online computation load of 10.9 s (no parallelization used). Similarly, for the case comprised of only the GPS satellites with elevation > 20°, 13 disjoint ones were estimated using our proposed ZSM with an online computation load of 5.5 s. We also observed that, in both cases, our proposed ZSM successfully detected the mode that contained the true receiver location with a low centroid error of [0.9,6.0] m and a bound of [21.4,55.1] m in local map coordinates for the 14 GPS satellite case. On the other hand, a higher centroid error of [1.6,10.8] m and bound of [21.4,122.5] m in local map coordinates was observed for the six GPS satellite case. Unlike in Section 4.4.3 where increasing the number of satellites did not significantly affect the characteristics of the mode containing the true receiver location, in this case, we observed an improvement in our ZSM’s performance. As explained before, this can be attributed to the further addition of useful geometric constraints in Figure 14(c) as compared to Figure 14(d) using the GPS shadows extracted from the remaining eight satellites.
5 CONCLUSIONS
We presented zonotope shadow matching (ZSM), a novel approach to set-valued position estimation for 3D-map-aided GNSS. The method is achieved by leveraging constrained zonotopes, a recent advance in polytope representation. We computed the 3D buildings and 2D GNSS shadows using constrained zonotopes, and then refined the coarse set-valued receiver position estimate based on if the receiver lay inside or outside the GNSS shadow (which is judged by C / N0 values). Using simulated experiments on a simple 3D building map and a dense 3D map of San Francisco, we validated that the proposed ZSM achieves a high point-valued estimation accuracy with high certainty (small position bounds). Importantly, while achieving a comparable point-valued accuracy as that of the conventional shadow matching (SM) technique, ZSM can compute a set-valued state estimate that is independent of the point-valued discretization. We also demonstrated the easy scalability of ZSM with the number of buildings considered. This indicates that ZSM is a promising method for providing robustness guarantees for safety-critical GNSS applications, since such set-valued estimates can be used as measurements for set-valued robust estimators.
HOW TO CITE THIS ARTICLE
Gao, G., Bhamidipati, S., & Kousik, S. (2022). Set-valued shadow matching using zonotopes for 3D-map-aided GNSS localization. NAVIGATION, 69(4). https://doi.org/10.33012/navi.547
ACKNOWLEDGMENTS
We want to thank the National Science Foundation (NSF) and the Air Force Research Laboratory (AFRL) for funding this work. We also want to thank Adyasha Mohanty and Asta Wu for reviewing our paper.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.