# Mosaic Zonotope Shadow Matching for Risk-Aware Autonomous Localization in
Harsh Urban Environments

Daniel Neamati, Sriramya Bhamidipati and Grace Gao  Stanford University

###### Abstract

Risk-aware urban localization with the Global Navigation Satellite System
(GNSS) remains an unsolved problem with frequent misdetection of the user’s
street or side of the street. Significant advances in 3D map-aided GNSS use
grid-based GNSS shadow matching alongside AI-driven line-of-sight (LOS)
classifiers and server-based processing to improve localization accuracy,
especially in the cross-street direction. Our prior work introduces a new
paradigm for shadow matching that proposes set-valued localization with
computationally efficient zonotope set representations. While existing
literature improved accuracy and efficiency, the current state of shadow
matching theory does not address the needs of risk-aware autonomous systems.
We extend our prior work to propose Mosaic Zonotope Shadow Matching (MZSM)
that employs a classifier-agnostic polytope mosaic architecture to provide
risk-awareness and certifiable guarantees on urban positioning. We formulate a
recursively expanding binary tree that refines an initial location estimate
with set operations into smaller polytopes. Together, the smaller polytopes
form a mosaic. We weight the tree branches with the probability that the user
is in line of sight of the satellite and expand the tree with each new
satellite observation. Our method yields an exact shadow matching distribution
from which we guarantee uncertainty bounds on the user localization. We
perform high-fidelity simulations using a 3D building map of San Francisco to
validate our algorithm’s risk-aware improvements. We demonstrate that MZSM
provides certifiable guarantees across varied data-driven LOS classifier
accuracies and yields a more precise understanding of the uncertainty over
existing methods. We validate that our tree-based construction is efficient
and tractable, computing a mosaic from 14 satellites in 0.63 seconds and
growing quadratically in the satellite number.

###### keywords:

GNSS, Shadow Matching, Set-based, Risk-Aware Localization

##  1 Introduction

A growing number of autonomous transportation systems are being developed for
safety-critical urban applications [1], including traffic management, lane
keeping, ride-share services, cargo transport, and medical-aid delivery. With
recent advances in sensor technologies and computing resources, there is an
increasing interest in using Artificial Intelligence (AI)-driven tools for
localization of autonomous vehicles, such as self-driving cars and unmanned
aerial vehicles (UAVs) [2, 3, 4]. Specifically, AI tools are popular in
localization frameworks that rely on data-driven classifiers to identify
distinct features in sensor measurements. For example, AI tools have been
successful in object detection and classification for computer vision [5] and,
similarly, for satellite LOS or NLOS classification for Global Navigation
Satellite System (GNSS) [6, 7]. However, a lack of widespread risk-aware and
certifiable localization inhibits the AI decision-making agents (in this case
autonomous systems) from making optimal decisions in safety-critical settings,
thus leading to potentially catastrophic consequences [8].

Of many possible sensor choices, GNSS receivers have emerged as a standard
choice for the outdoor localization of many autonomous systems since these
receivers provide absolute user location [9]. However, many tall buildings in
close proximity form an urban canyon that degrades the performance of GNSS-
based localization, often incorrectly locating the user’s street or side of
the street. Specifically, tall buildings block, reflect, scatter, and diffract
the direct Line-of-Sight (LOS) GNSS signals, thus inducing Non-Line-of-Sight
(NLOS) and multipath effects [10]. The broad field of 3D Map-Aided GNSS (3DMA-
GNSS) uses 3D building maps and LOS classifiers to detect, mitigate, and
correct the NLOS and multipath effects that degrade GNSS performance [11, 12,
13]. A 3D building map encodes the 3D structure of the city’s infrastructure
and a LOS classifier analyzes signal properties to predict the probability
with which the user is in line of sight of a satellite. However, the
complexity of the real-world impacts both the realism of 3D building maps and
the accuracy of LOS/NLOS classifiers. For example, differences in the material
of the building’s facade, which are difficult to include in the 3D building
map, change the way that the GNSS signals reflect, scatter, and diffract off
the building [14, 15]. The local environment also changes across seasons and
with dynamic obstacles. Various methods augment 3D building maps with digital
surface models, antenna hardware advances, on-board imaging, satellite
imagery, and LiDAR [16, 17, 18, 19]. However, current 3DMA-GNSS works struggle
to meet the computational efficiency needed to include these real-time 3D
building map augmentations on board. New data-driven/AI-driven LOS classifiers
employ decision trees, support vector machines, or neural networks on many
signal features including carrier-to-noise density ratio ($C/N_{0}$) and
pseudorange residuals to estimate the probability that the user is in line of
sight of the satellite [6, 20]. However, environmental complexities make
classification difficult, with past data-driven works reporting accuracies of
roughly 80% to 90% [6, 20].

Advancements in computational capabilities and increasing numbers of GNSS
signals have steadily helped improve both the 3D building map processing and
AI-driven LOS classifier accuracy. The resources for computation processing
(e.g., RAM) and storage memory (e.g., flash memory) have increased
considerably since 3DMA-GNSS was first proposed. In 2013, authors used
previous top-of-the-line smartphones, such as the Samsung Galaxy S3 with up to
2 GB RAM and 16 GB to 64 GB flash memory [21]. In comparison, the recent
Samsung Galaxy S22 Ultra has up to 12 GB RAM and 128 GB to 1 TB flash memory.
On-board processing on autonomous systems, such as cars and drones, has
correspondingly improved. In 2013, authors from [21, 22] primarily used GPS
and GLONASS signals from one frequency band (e.g., GPS L1 C/A signal). More
modern antennas in phones and autonomous systems now use multiple recently-
available constellations, such as BeiDou and Galileo, and multiple
frequencies, such as the GPS L1 and L5 signals [23].

At the methodological level, there are two primary categories of 3DMA-GNSS:
ray-tracing and shadow matching. That said, we also acknowledge that recent
end-to-end machine learning (ML) methods from Google [24] take a fundamentally
different approach to 3DMA-GNSS, though this approach has not exhibited much
available literature and we cannot assess its risk-awareness. The first
primary category is ray-tracing, which considers a large sample of possible
signal paths between the user and the satellite. By computing possible
scattering and reflections off buildings, ray-tracing uses the computed
distribution in signal path length to recover a corrected GNSS signal. For a
sufficient number of samples, ray-tracing provides an accurate and precise
understanding of the distribution of possible signals and the associated
uncertainties [12, 25, 26]. However, ray-tracing is too computationally
expensive even for a single satellite and cannot be used in real-time
applications without incurring large server overhead.

The second primary approach to 3DMA-GNSS is shadow matching. Grid-Based Shadow
Matching (GBSM) is the current state-of-the-art for shadow matching and is
often mixed with 3DMA GNSS ranging in the Intelligent Urban Positioning (IUP)
methodology [27, 28, 29]. GBSM starts with a grid overlaid over the entire
city at a fixed resolution where each grid point represents a possible
receiver location. GBSM derives an offline building boundary look-up table for
each grid point that maps elevation and azimuth with prespecified angular
discretization to satellite visibility. Specifically, a satellite would be
visible at a given elevation and azimuth if no buildings in the 3D building
map would block that satellite’s LOS signal. The look-up table method provides
online efficiency for a fixed map across large numbers of satellites and has
been consistently demonstrated in real-time settings [13, 21, 27, 28, 29, 30,
31]. Unlike ray-tracing, shadow matching can provide localization even when
the user does not receive a GNSS signal from a satellite (e.g., if the
buildings fully blocked the satellite). Works such as [6, 20] have
successfully incorporated AI-driven classifiers with the shadow matching
building boundary look-up tables to improve localization. Hence, shadow
matching has been a successful method to improve urban GNSS accuracy in real-
time.

However, urban GNSS localization accuracy is not the only objective. Safety-
critical autonomous systems require further guarantees on the urban
localization solution to ensure the understanding of risk necessary for
downstream AI decision-making [8]. From a bigger picture outlook, this gap is
not unique to GNSS localization. Authors in [8] present a thorough survey of
uncertainty quantification in AI and learning-based tools. A large volume of
work has addressed uncertainty for single specific components such as
uncertainty quantification in computer vision or reinforcement learning.
However, we must also provide the theoretical tools for transitioning
uncertainty quantification in perception systems to AI decision-making agents
in a risk-aware manner. As throughout risk-aware AI, the push toward risk-
aware urban GNSS localization requires new theoretical developments. To reach
risk-aware 3DMA-GNSS, one must meet the following two objectives: (1)
formulate precise models of the uncertainty in real-time while tractably
scaling with increasing numbers of available signals and (2) provide
guarantees on the uncertainty bounds of the urban localization solution across
a vast space of AI-driven LOS classifier designs.

Several works in GBSM have considered the implications of these risk-aware
objectives on 3DMA-GNSS. In [21, 22], the authors discussed their method to
compute the building boundaries, which can compute a 1-meter grid resolution
for a 500 m by 500 m AOI in less than 4 days. This method of computing
building boundaries cannot adapt to real-time changes. Moreover, with the
grid-based approach, one cannot assess if two highly weighted grid points use
the same or different satellite visibility information. This raises the
concern that a fixed number of top-scoring output grid cells could contain
identical misclassification, thereby conflating the shadow matching
uncertainty with the grid resolution. GBSM also did not provide self-
monitoring on the initial localization. Instead, they assumed large areas of
interest to begin processing and temporally filter from the initial
localization assuming a lock on the user location. While many recent works in
IUP use a grid-based filter to combine information and handle the highly non-
Gaussian uncertainties in urban localization, many of the output statistics
use a weighted average approach with an implicit Gaussian approximation.
Lastly, the fixed grid spacing and fixed elevation/azimuth spacing enforce a
fixed precision of the overall GBSM approach. Improving the precision requires
finer grid spacing and finer elevation/azimuth spacing which come at an
increased computational cost. To the authors’ knowledge, GBSM cannot yet
provide certifiable guarantees on the user uncertainty bounds. In our results
(Section 3), we illustrate the how the lack of risk-awareness in GBSM can
impede an accuracy and precision understanding of the shadow matching
uncertainty.

Figure 1: Our prior work on ZSM introduced set-based shadow matching, which
efficiently extracts a GNSS shadow (magenta) using set-valued representations,
namely constrained zonotopes [32]. We denote the building in gray, the shadow
direction with the black thick arrow, the shadow volume in cyan, the ground
plane in brown, and the extracted GNSS shadow in magenta. ZSM recursively
refines the area of interest with GNSS shadows, but assumes an ideal LOS
classifier, and is thus brittle against any LOS classifier inaccuracies.
(Figure adapted from [32]).

Our recent work proposes a new shadow matching paradigm: set-based shadow
matching [32]. A set-valued approach computes a continuum of state estimates
rather than considering a grid of position candidates. In related set-based
works, authors use an initial set of states and measurement bounds to provide
certifiable guarantees on the state estimation by checking if the final
estimated set lies within user-specified safety bounds [33, 34]. In literature
[33, 35, 36], many set-valued representations exist for parameterizing the
states. The zonotope is a widely popular set-representation in robotics for
computing reachable sets of dynamical systems that enable formally-verified
path planning and collision avoidance [35]. Zonotopes are convex, symmetrical
polytopes that can propagate set-valued state estimates using fast vector
concatenation operations. Furthermore, constrained zonotopes generalize
zonotopes to represent any convex polytopes, thereby avoiding the limitation
of symmetry [34]. Our prior work [32] developed a novel Zonotope Shadow
Matching (ZSM) algorithm that computes the set-based GNSS shadow based on the
3D building map, as shown in Figure 1. ZSM parameterizes each
satellite/building pair with constrained zonotopes, extends the building along
shadow direction from the GNSS satellite, and performs set intersection with
the ground plane. Thereafter, this ZSM work iteratively refines an area of
interest (AOI) with an ideal LOS classifier to judge if the receiver is inside
or outside the GNSS shadow; if outside, set subtraction is executed, otherwise
set intersection. While this work demonstrates great success in achieving
positioning accuracy with high certainty for various simulation scenarios,
there still exists a major shortcoming: ZSM is brittle to misclassification.
In iteratively processing the GNSS shadows, even a single incorrect LOS/NLOS
classification can lead to divergence of the algorithm, resulting in incorrect
localization, i.e., a null set for the user position. Simply, ZSM alone is not
risk-aware.

In this work, we extend the set-based shadow matching paradigm to allow for
risk-aware 3DMA-GNSS. Specifically, we propose Mosaic Zonotope Shadow Matching
(MZSM), which uses the shadow extraction portion of ZSM as a front-end
processing unit and a novel polytope mosaic as a back-end risk-aware
localization unit. Effectively, we use the portion of ZSM that extracts set-
based GNSS shadows from the GNSS satellite position and 3D building map. Our
novel polytope mosaic stitches the GNSS shadows and their complements (the
rest of the AOI not within GNSS shadow) with the LOS probabilities from the
AI-driven LOS classifier into a risk-aware framework. MZSM exhaustively
considers all possible classifier misclassifications to guarantee uncertainty
bounds on the output set-based localization solution. To ensure that this
exhaustive mapping of classifier misclassifications is tractable and
efficient, we employ a memory-efficient binary tree that grows quadratically
with the number of satellites rather than the exponential growth of fully-
branching binary trees. To ensure MZSM can operate with off-the-shelf AI-
driven classifiers and adapt to future classifier designs, we take a novel
classifier-agnostic approach. Lastly, we show how our tree-based architecture
can predict the expected value shadow matching distribution based on the
classifier accuracy statistics. This prediction capability will allow AI
decision-making agents to consider the expected uncertainty on the shadow
matching distribution in advance, which enables risk-awareness in downstream
autonomous planning. Altogether, we aim to further the state-of-the-art in
shadow matching to include risk-aware urban localization via MZSM.

###  1.1 Contributions and Paper Organization

Our key contributions are as follows:

  1. 1.

We develop the theoretical foundations for a polytope mosaic that sequential
processes LOS classifier probabilities and GNSS shadow polytopes for each GNSS
satellite to provide classifier-agnostic certifiable guarantees on the user
position.

  2. 2.

We design a recursively expanding binary tree that performs set operations to
bifurcate the area of interest into smaller polytopes at each tree level, and
thereafter, probabilistically weighs each tree branch based on the LOS
classifier probability. To ensure tractability, we special case different
shadow overlap conditions to produce a memory-efficient binary tree. Our
method computes an exact shadow matching distribution from which we determine
a confidence collection of smaller polytopes that collectively satisfies a
user confidence level.

  3. 3.

We produce an analysis method to estimate the expected shadow matching
distribution from a classifier’s accuracy statistics (more broadly, its
uncertainty quantification statistics). We demonstrate that this analysis
method arises from the information organization in our binary tree and the
set-based localization in our mosaic. We discuss how a system designer and an
AI decision-making agent can leverage our analysis method to forecast the
shadow matching mosaic in design or online planning.

  4. 4.

We experimental validate our risk-aware shadow matching focusing on the
improvements of classifier-agnostic evaluation, set-based bounds, and
efficiency with tree-based processing.

The remainder of the paper is organized as follows. Section 2 describes our
proposed MZSM implementation through the tree-based MZSM construction. We end
Section 2 with a discussion of set-based uncertainty bounds and an
illustrative example. We validate our key risk-aware contributions of MZSM in
Section 3 via simulated experimental results. Finally, Section 4 provides
concluding remarks.

##  2 Proposed Mosaic Zonotope Shadow Matching Algorithm

We first describe our model and assumptions of the GNSS satellite signals, the
initial area of interest, and the 3D building map. We then summarize the
procedure for extracting GNSS shadows, which is based on our prior work [32].
We finally explain the details of our tree data structure that takes in the
GNSS shadows as input to construct the exact shadow matching distribution
(i.e., no discretization), which provides certifiable guarantees on the
receiver position. Figure 2 provides a roadmap for the methods.

Figure 2: Our MZSM algorithm with two key modules: a) extraction of GNSS
shadows, which is based on our prior work [32]; and b) construction of a
polytope mosaic. We take the 3D building map, the area of interest, and the
known GNSS satellite positions to compute a set collection of GNSS shadow
polytopes. We construct our polytope mosaic by taking in the area of interest,
the LOS classifier probabilities, and the set of GNSS shadow polytopes. Our
MZSM provides risk-aware 3DMA-GNSS localization, wherein the outputs include a
confidence collection set of receiver positions that collectively satisfy a
user confidence level and a PDF of user position over the mosaic conditioned
on the AOI.

###  2.1 Notation

Throughout, scalars are italic lowercase. Sets and sets of sets are in script
uppercase. We use a subscript for indexing and a text superscript for direct
labeling. Probability is denoted $p(\cdot)$. We use the shorthand nation
$\\{\cdot_{i}\\}_{i=1}^{n}$ to construct a set of $n$ elements where each
element is indexed with $i$.

###  2.2 System Inputs

MZSM takes in three inputs: (1) an AOI, (2) GNSS shadow polytopes, and (3) LOS
probabilities from the classifier. We briefly discuss each input below.

First, the AOI is an initial (often coarse) bounded estimate of the 2D
receiver location. As in past works [32, 37], we use a terrain model from the
3D building map to enforce that the 2D AOI properly accounts for the city
topography. Mathematically, we parameterize the AOI ($\mathcal{A}$) by either
a single 2D polytope or a collection of 2D polytopes to match the terrain.

$\mathcal{A}=\\{A_{k}\\}_{k=1}^{n_{\mathrm{\textnormal{\tiny{AOI}}}}}\subseteq\mathcal{G},$
(1)

where $A_{k}$ denotes a single polytope among
$n_{\mathrm{\textnormal{\tiny{AOI}}}}$ total that form $\mathcal{A}$, and
$\mathcal{G}$ represents the 2D ground plane. We remove the building
footprints from the AOI, effectively leaving the polytopes with holes. If no
prior estimate is available, we can set the AOI ($\mathcal{A}$) equal to the
ground plane ($\mathcal{G}$), where $\mathcal{G}$ is considered to be bounded.
Since the AOI is a coarse estimate, risk-aware shadow matching must self-
monitor the AOI to determine if the AOI is consistent with the incoming GNSS
information. Our MZSM method is the first to incorporate this AOI self-
monitoring feature (elaborated upon in Section 2.4 and quantitatively assessed
in Section 3.6).

Second, we compute the GNSS shadows using the 3D building map per the method
in our prior ZSM work [32]. We briefly summarize the important aspects of ZSM
for MZSM. Standard 3D city maps usually comprise a collection of 3D vertices.
We convert the point-based 3D city map to a constrained zonotope city map via
triangulation, as described further in [32]. We assume that the user has
access to measurements from $n_{\mathrm{\textnormal{sat}}}$ GNSS satellites.
We denote the observed set of GNSS satellites by
$\mathcal{S}=\\{s_{j}\\}_{j=1}^{n_{\mathrm{\textnormal{sat}}}}$, where each
$s_{j}\in\mathbb{R}^{3}$ represents the satellite location and
$n_{\mathrm{\textnormal{sat}}}$ is the total number. For each satellite-
building pair, we compute the shadow volume (cyan in Figure 1) and project the
shadow onto the ground plane. We convert the constrained zonotope shadow for
each satellite-building pair into a vertex representation and take the union
across buildings. This yields one GNSS shadow ($\mathcal{C}_{j}$, magenta in
Figure 1) as a vertex representation for each satellite $j$. Full details on
the shadow extraction algorithm and online processing capabilities are
available in [32]. Importantly for the mosaic, the GNSS shadows are completely
independent of the receiver location. In other words, the shadows are fixed
with respect to the buildings and satellites such that they are the same for
all receivers in the scene and change only as the satellite orbits.

Third, the user must determine the probability with which they are in the LOS
of the satellite or in the satellite’s GNSS shadow. We leverage LOS/NLOS
classifiers to determine the probability $p(s_{j}=\mathrm{\textnormal{LOS}})$
that the user is in the LOS of the $j$-th satellite, $\forall
s_{j}\in\mathcal{S}$. By extension, we express the probability that the user
is in the 2D GNSS shadow as $p(s_{j}=\mathrm{\textnormal{NLOS}})$, where
$p(s_{j}=\mathrm{\textnormal{LOS}})=1-p(s_{j}=\mathrm{\textnormal{NLOS}})$.
Our proposed algorithm is independent of the choice of classifier utilized,
but, in practice, system designer can largely select any classifier from
literature [6, 20].

###  2.3 Tree-based processing and expansion from a collection of GNSS shadows

We now proceed to build a risk-aware polytope mosaic from the AOI denoted by
$\mathcal{A}$, the set of GNSS shadows denoted by
$\\{\mathcal{C}_{j}\\}_{j=1}^{n_{\mathrm{\textnormal{sat}}}}$, and the LOS
classifier probability values denoted by
$p(s_{j}=\mathrm{\textnormal{LOS}}),\,\forall s_{j}\in\mathcal{S}$. We
sequentially process the extracted GNSS shadow polytopes $\mathcal{C}_{j}$ to
split the AOI polytope $\mathcal{A}$ into various smaller polytopes. Together
the smaller polytopes form a mosaic over the AOI. Simply as a polytope mosaic,
it is not straightforward to estimate which of the smaller polytopes is the
receiver most likely to be in. To address this, we design a probabilistic
polytope mosaic that sequentially processes the extracted GNSS shadow
polytopes $\mathcal{C}_{j}$ with the LOS/NLOS probabilities
$p(s_{j}=\mathrm{\textnormal{LOS}})$ or $p(s_{j}=\mathrm{\textnormal{NLOS}})$
associated with the $j$-th GNSS satellites. Thus, each smaller polytope has a
unique identification based on the associated joint LOS/NLOS probabilities
from the set of all observed GNSS satellites $\mathcal{S}$.

Figure 3: Geometry of the three overlap cases: no shadow overlap, split
overlap, and full shadow overlap. For each case, we illustrate how the overlap
between the shadow (blue) with the parent node polytope (green) dictates the
children node polytopes. The no shadow overlap case generates an LOS child
polytope (orange) but no NLOS child polytope (purple). Equivalently, the NLOS
child polytope is the empty set. The full shadow overlap case generates an
NLOS child but no LOS child. Equivalently, the LOS child polytope is the empty
set. Only the split overlap case generates non-empty LOS and NLOS children.
The split overlap can generate single polytopes (such as the either of the
NLOS child polytope), disjoint sets (such as the left LOS child polytope), or
sets with holes (such as the right LOS child polytope).

We organize these sequential operations as a recursively-expanding binary tree
to connect information across satellites and ensure computational efficiency.
The tree root is the AOI polytope $\mathcal{A}$ with probability
$p(\mathcal{A})$ that the user is in the AOI. Generally, $p(\mathcal{A})=1$.
We discuss how we monitor if the AOI is violated later in the methods (Section
2.4). Each node in the tree encodes a set-based location and associated
probability that the user is in that set. Each set-based location is a
polytope or a collection of disjoint polytopes, which is collectively denoted
by $\lozenge$. The associated probability is denoted by $p(\lozenge)$. Each
tree branch represents a LOS- or an NLOS-based bifurcation of the parent node
polytope into two children nodes. The LOS/NLOS bifurcations and the root AOI
uniquely describe a node in the tree. A node with two empty subtrees (no
further branching) is called a leaf. Each new GNSS satellite allows us to
expand the tree, such that the maximum depth of the tree equals the number of
GNSS satellites $n_{\mathrm{\textnormal{sat}}}$.

When incorporating another GNSS shadow to expand a depth $(j-1)$ tree to the
subsequent depth $j$ tree, there are three possible overlap cases between the
GNSS shadow and a node polytope. These three overlap cases illustrated in
Figure 3 are listed as follows: (1) the GNSS shadow and parent node polytope
do not overlap, (2) the GNSS shadow either partially or fully lies within the
parent node polytope (thereby causing a split in the parent node polytope),
and (3) the parent node polytope fully lies within the GNSS shadow. We can
express the overlap conditions mathematically via Eqs. (2)-(4) where we label
$\lozenge_{j-1,n}$ as the $n$-th node polytope at the $(j-1)$-th tree level.
The $0$-th node polytope at the $0$-th tree level represents the tree’s root,
and thus $\lozenge_{0,0}$ equals $\mathcal{A}$.

No Shadow Overlap:
$\displaystyle\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\emptyset$ (2) Split
Overlap:
$\displaystyle(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\emptyset)\text{ and
}(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\lozenge_{j-1,n})$ (3) Full Shadow
Overlap: $\displaystyle\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\lozenge_{j-1,n}$
(4)

Note that the condition specified for full shadow overlap in Eq. (4)
essentially implies $\lozenge_{j-1,n}\subseteq\mathcal{C}_{j}$. These three
cases provide a distinct and complete (though non-unique) categorization of
all possible overlap cases (Proposition 1).

######  Proposition 1.

In branching any $(j-1)$-th tree level to the next based on the $j$-th GNSS
shadow $\mathcal{C}_{j}\neq\emptyset$, the outcome of overlaying
$\mathcal{C}_{j}$ onto any parent node polytope
$\lozenge_{j-1,n}\neq\emptyset,~{}\forall n,$ in the $(j-1)$-th tree level can
be uniquely characterized as either (1) no shadow overlap, (2) split overlap,
or (3) full shadow overlap. That is, these cases are distinct and complete.

###### Proof.

See A. ∎

As a result of the overlap between a GNSS shadow and any position node
polytope designed as a parent, we generate a maximum of two children node
polytopes in the next tree level. One child, designated as LOS child,
represents the portion of the parent node polytope that is in the LOS of the
satellite. The other child, namely NLOS child, indicates the portion that is
not in the LOS of the satellite or equivalently, the portion that is in the
satellite’s GNSS shadow. As in [32], we calculate the LOS case with set
subtraction in Eq. (5) and the NLOS case with set intersection in Eq. (6).
However, the existence (not empty) or feasibility (non-zero probability) of
these children node polytopes depend on the overlap condition (from
Proposition 1) satisfied by the parent node and the GNSS shadow.

$\displaystyle\lozenge^{\mathrm{\textnormal{LOS child}}}_{j,m}$
$\displaystyle=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}-\mathcal{C}_{j}=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\cap\mathcal{C}_{j}^{\mathsf{c}}$
(5) $\displaystyle\lozenge^{\mathrm{\textnormal{NLOS child}}}_{j,l}$
$\displaystyle=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\cap\mathcal{C}_{j}$
(6)

where $\lozenge^{\mathrm{\textnormal{LOS child}}}_{j,m}$,
$\lozenge^{\mathrm{\textnormal{NLOS child}}}_{j,l}$, and
$\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}$ denotes the node polytopes,
such that $\lozenge^{\mathrm{\textnormal{LOS child}}}_{j,m}$,
$\lozenge^{\mathrm{\textnormal{NLOS child}}}_{j,l}$ are the children nodes of
the parent node $\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}$.
Additionally, $\mathcal{C}_{j}^{\mathsf{c}}$ represents the complement of the
$j$-th GNSS shadow polytope with
$\mathcal{C}_{j}^{\mathsf{c}}=\mathcal{A}-\mathcal{C}_{j}$. Eqs. (5) and (6)
recursively produce the children from the parent.

The no shadow overlap and full shadow overlap cases generate empty child
nodes. When the GNSS shadow and the parent node polytope do not overlap,
$\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\cap\mathcal{C}_{j}=\emptyset$,
as seen in Eq. (2) and the no shadow overlap case of Figure 3. Thus
$\lozenge^{\mathrm{\textnormal{LOS
child}}}_{j,m}=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}$ and
$\lozenge^{\mathrm{\textnormal{NLOS child}}}_{j,l}=\emptyset$. On the other
hand, when the parent node polytope fully overlaps with the GNSS shadow
polytope,
$\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\cap\mathcal{C}_{j}=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}$,
as seen in Eq. 4 and the full shadow overlap case of Figure 3. Thus,
$\lozenge^{\mathrm{\textnormal{LOS child}}}_{j,m}=\emptyset$ and
$\lozenge^{\mathrm{\textnormal{NLOS
child}}}_{j,l}=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}$. During a
split overlap case seen in Eq. (3), both children node polytopes exist as long
as $\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\neq\emptyset$ and
$\mathcal{C}_{j}\neq\emptyset$. Intuitively, the empty children nodes in any
of these two cases imply redundant or inconsistent measurements. For instance,
consider the no shadow overlap case at any $j$-th tree level. Given that the
receiver lies within the parent node polytope, no overlap means that all
locations within the parent node polytope are in LOS of the $j$-th GNSS
satellite. The LOS child node thereby reflects redundant information while the
NLOS child reflects incompatible information. The reverse relationship holds
true for full shadow overlap, i.e., the LOS child node reflects incompatible
information, and the NLOS child node reflects redundant information.

We assign a probability score to the children corresponding to
$p(s_{j}=\mathrm{\textnormal{LOS}})$ for the LOS child node polytope and
$p(s_{j}=\mathrm{\textnormal{NLOS}})=1-p(s_{j}=\mathrm{\textnormal{LOS}})$ for
the NLOS child node polytope. Starting with a pre-defined value of
$p(\lozenge^{\mathrm{\textnormal{parent}}}_{0,0})=p(\mathcal{A})$, we arrive
at recursive equations (Eq. (7)-(8)) with conditional independence between the
GNSS satellite signals. The probability scores decrease monotonically for
nodes deeper in the tree, since probabilities are between zero and one.

$\displaystyle p\Big{(}\lozenge^{\mathrm{\textnormal{LOS
child}}}_{j,m}\Big{)}$ $\displaystyle=p(s_{j}=\mathrm{\textnormal{LOS}})\cdot
p\Big{(}\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\Big{)}$ (7)
$\displaystyle p\Big{(}\lozenge^{\mathrm{\textnormal{NLOS
child}}}_{j,l}\Big{)}$ $\displaystyle=p(s_{j}=\mathrm{\textnormal{NLOS}})\cdot
p\Big{(}\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\Big{)}$ (8)

Starting with the AOI polytope, we sequentially process the GNSS shadows and
recursively expand the tree to arrive at a completed binary tree. The leaves
of this tree form a probabilistic mosaic of polytopes. Our binary tree is
commutative, the mathematical proof of which is discussed in Proposition 2.
Intuitively, commutativity implies that we arrive at the same completed binary
tree irrespective of the order in which we overlay the GNSS shadows onto the
AOI. This is a crucial property in ensuring risk-aware 3DMA-GNSS localization
since an observation error from the AI-driven LOS classifier in the LOS/NLOS
probability of a GNSS shadow processed early in the ordering will not disrupt
the tree at later operations. That is, we preserve the conditional
independence between GNSS signals.

######  Proposition 2.

Given a set of GNSS shadow polytopes
$\\{\mathcal{C}_{j}\\}_{j=1}^{n_{\mathrm{\textnormal{sat}}}}$ and an AOI
polytope $\mathcal{A}$, the leaf polytopes (node polytopes with no further
branching) generated from our completed binary tree are the same for any
ordering of the AOI and the GNSS shadows.

###### Proof.

See B. ∎

###  2.4 Memory-efficient binary tree for a tractable mosaic

During the no shadow overlap and full shadow overlap cases explained earlier,
the NLOS child and the LOS child polytope are empty, respectively. Expanding
the tree with an empty set as the parent node polytope leads to empty sets as
the children nodes as well. Mathematically, if
$\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}=\emptyset$, then
$\lozenge^{\mathrm{\textnormal{LOS
child}}}_{j+1,l}=\emptyset-\mathcal{C}_{k}=\emptyset$ via Eq. (5) and
$\lozenge^{\mathrm{\textnormal{NLOS
child}}}_{j+1,m}=\emptyset\cap\mathcal{C}_{k}=\emptyset$ via Eq. (6). In
Proposition 3 we demonstrate that the non-empty leaf polytopes are sufficient
to produce the complete mosaic over the AOI.

######  Proposition 3.

Given a set collection of leaf polytopes, the union of the non-empty leaf
polytopes exactly matches the tree’s root, i.e., AOI polytope ($\mathcal{A}$).

###### Proof.

See C. ∎

However, the probabilities of the empty parent and empty children are non-
zero. These probabilities represent the possibility that one of the parent’s
branches is not compatible. The probability is sizeable when the AOI does not
include the true receiver location. Thus, we do not need to propagate an empty
set in the tree from a geometric sense, but we do need to keep track of the
total probability of the empty nodes to self-monitor the probability that our
AOI is violated. We use this to our advantage: we avoid expanding the tree
when we detect that a non-empty parent node will yield an empty child
polytope. Instead of tree expansion, we update the probability score of the
parent node with on LOS/NLOS probability of the non-empty child node polytope.
In particular, we identify the empty children nodes under each overlap
condition by assessing the area estimates of the LOS child and NLOS child. We
then update the probabilities using the same method as Eqs. (7) and (8), but
assigning the probability update to the parent. We demonstrate this self-
monitoring capability across classifiers later in Section 3.6 (Figure 11).

We primarily exclude the empty sets to maintain a low and tractable memory
footprint for the tree. For instance, if we divided every parent node polytope
into two children nodes (without excluding empty sets) for each GNSS shadow
polytope, the binary tree would grow to have $2^{n}$ leaf polytopes whose sum
of probability scores satisfy
$\sum\limits_{n}p(\lozenge_{n_{\mathrm{\textnormal{sat}}},n})=\sum\limits_{n}p(\lozenge^{\mathrm{\textnormal{leaf}}}_{n})=p(\mathcal{A})$.
Geometrically, the total number of leaf polytopes will be far fewer than the
$2^{n}$ size for a full binary tree if we exclude the empty sets. This
mathematically follows from the Lazy Caterer’s sequence (integer sequence
A000124 of OEIS [38]) where the number of polygonal sections formed from $n$
cuts grows quadratically with the number of cuts. Intuitively, given a
sufficiently deep tree, processing the next GNSS shadow polytope cannot result
in the split overlap case for all the leaf nodes. Instead, the GNSS shadow
will fully overlap many nodes and will not overlap with many other nodes. We
experimentally validate this quadratic complexity later in Section 3.6
(Figures 12 and 13).

Lastly, we must ensure that our memory-efficient binary tree maintains the
risk-aware self-monitoring of the full binary tree. In our memory-efficient
binary tree, the AOI violation probability (i.e., probability of inconsistent
information) is directly calculable from the root and leaf probabilities with
Eq. (9).

$p_{\emptyset}=p(\mathcal{A})-\sum\limits_{l}p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l})$
(9)

where $\sum\limits_{l}p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l})$ denotes
the sum of probability scores at the non-empty leaf nodes in the memory-
efficient binary tree. By construction, $p_{\emptyset}$ is the total
probability of all the empty leaf nodes in the full tree. We can use
$p_{\emptyset}$ to diagnose the posterior confidence in our AOI. If
$p_{\emptyset}\approx 1-p(\mathcal{A})$, then the mismatch in $p_{\emptyset}$
is more attributed to inaccuracies in the LOS classifier. For example, given
$4$ GNSS satellite signals with
$p(s_{j}=\mathrm{\textnormal{NLOS}})=0.1~{}\forall j\in\\{1,\cdots,4\\}$, the
leaf polytope that represents lying within the GNSS shadow polytope four
consecutive times has a probability score of $(0.1)^{4}=0.0001$. If this
position leaf polytope is an empty set (i.e., area equals zero), then either
the AOI is incorrect or at least one of the satellites should have a higher
LOS probability. Further, if this leaf with four NLOS branches is the only
empty leaf polytope in the entire binary tree, then $p_{\emptyset}=0.0001$. On
the other hand, if $p_{\emptyset}\approx p(\mathcal{A})$, then the mismatch is
more attributed to an incorrect AOI. To understand this, let’s consider a
different case for the same example as earlier with $4$ GNSS satellites.
Suppose that the position leaf polytope that represents lying within the
complement of the GNSS shadow polytope four consecutive times is an empty set.
However, the probability score of this position leaf polytope exhibits a
significant value of $(0.9)^{4}=0.6561$. In this case, we are more likely to
trust that the GNSS signals were classified correctly as LOS, but the initial
AOI polytope did not include the true receiver position. If this leaf with
four LOS branches is the only empty leaf polytope in the entire binary tree,
then $p_{\emptyset}=0.6561$. Hence, at low $p_{\emptyset}$ values the tree’s
self-monitor is confident in the AOI. At moderate to high $p_{\emptyset}$
values the tree’s self-monitor is skeptical of the AOI and the AOI should be
enlarged accordingly. Therefore, our memory-efficient binary tree maintains a
worst-case quadratic complexity, only stores information of non-empty nodes,
and features built-in self-monitoring on the AOI root node.

We summarize the recursive expansion process to build the memory-efficient
binary tree in Algorithm 1.

1 Function  _recursive_tree_expansion(_node k,n, $\mathcal{C}_{j}$,
$p(s_{j}=\mathrm{\textnormal{LOS}})$_)_

2 /* Inputs: (1) nodek,n, (2) a GNSS shadow $\mathcal{C}_{j}$, and (3)
probability of LOS $p(s_{j}=\mathrm{\textnormal{LOS}})$ for GNSS satellite
$j$. The mosiac tree node stores: $\lozenge_{k,n}$, $p(\lozenge_{k,n})$, its
LOS_child node, and its NLOS_child node. */

3 if  _split overlap between $\lozenge_{k,n}$ and $\mathcal{C}_{j}$ (Eq. (3))_
then

4 if  _node k,n has no LOS_child and NLOS_child_ then

// Expanded tree at this node

5 add LOS_child with polytope per Eq. (5) and probability per Eq. (7)

6 add NLOS_child with polytope per Eq. (6) and probability per Eq. (8)

7 else

8 recursive_tree_expansion(_LOS_child of $\lozenge_{k,n},\ \mathcal{C}_{j},\
p(s_{j}=\mathrm{\textnormal{LOS}})$_)

9 recursive_tree_expansion(_NLOS_child of $\lozenge_{k,n},\ \mathcal{C}_{j},\
p(s_{j}=\mathrm{\textnormal{LOS}})$_)

10

11 else if  _no shadow overlap between $\lozenge_{j-1,n}$ and
$\mathcal{C}_{j}$ (Eq. (2))_ then

12 recursively propagate the LOS score ($p(s_{j}=\mathrm{\textnormal{LOS}})$)
to leaves (no further expansion)

13 else

// Remaining case is full shadow overlap (Eq. (4))

14 recursively propagate the NLOS score
($1-p(s_{j}=\mathrm{\textnormal{LOS}})$) to leaves (no further expansion)

return $\lozenge_{k,n}$ // updated mosaic tree node

15

Algorithm 1 Memory-Efficient Recursive Binary Tree Expansion (Snapshot)

###  2.5 MZSM Outputs: Shadow Matching Mosaic, AOI-conditioned PMF, and PDF
over the Mosaic

The most immediate output of the memory-efficient binary tree is the shadow
matching mosaic ($\mathcal{M}$, Eq. (10)), which denotes the set of all
$n_{\mathrm{\textnormal{leaves}}}$ leaf polytopes in the memory-efficient
binary tree.

$\mathcal{M}=\\{\lozenge^{\mathrm{\textnormal{leaf}}}_{l}\\}_{l=1}^{n_{\mathrm{\textnormal{leaves}}}}$
(10)

From $\mathcal{M}$, we produce an overall probability mass function (PMF) over
the mosaic ($\mathcal{M}$) and a separate self-monitor of AOI violation
($p_{\emptyset}$, Eq. (9)). Mathematically, we can condition the leaf node
probabilities $p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l})$ for each $l$-th
leaf to the AOI $p(\lozenge^{\mathrm{\textnormal{leaf}}}_{j}\mid\mathcal{A})$
as shown in Eq. (11). Together, the
$p(\lozenge^{\mathrm{\textnormal{leaf}}}_{j}\mid\mathcal{A})$ probabilities
from the PMF over $\mathcal{M}$.

$p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l}\mid\mathcal{A})=\frac{p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l})}{\sum_{l=1}^{n}p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l})}$
(11)

The user requires a probabilistic guarantee that they are within a bounded
region at a desired confidence to ensure risk-aware localization. From the
mosaic $\mathcal{M}$ and the PMF over the mosaic
($p(\lozenge^{\mathrm{\textnormal{leaf}}}_{j}\mid\mathcal{A})$), we can
generate two outputs that achieve this goal: (1) a probability density
function over the user location and (2) a set-based confidence collection. To
have certifiable guarantees, the latter is the preferred output.

In certain cases, we may need to work with the user location directly. In
position space, each leaf polytope is a uniform distribution
($\mathcal{U}_{\lozenge^{\mathrm{\textnormal{leaf}}}_{l}}$) with an irregular
shape. Namely, the position $(x,y)$ only has non-zero probability in the
polytope’s enclosed area, where $x$ and $y$ are any two orthogonal 2D
coordinates (e.g., North and East or cross-street and along-street). We
normalize the weights to produce a well-defined probability density function
(PDF) over the mosaic denoted by $f_{\mathcal{M}}(x,y)$, per Eq. (12).

$f_{\mathcal{M}}(x,y)=\left(\sum_{l=1}^{n_{\mathrm{\textnormal{leaves}}}}\mathcal{U}_{\lozenge^{\mathrm{\textnormal{leaf}}}_{l}}(x,y)\right)\left(\iint_{\text{AOI}}\sum_{l=1}^{n_{\mathrm{\textnormal{leaves}}}}\mathcal{U}_{\lozenge^{\mathrm{\textnormal{leaf}}}_{l}}(x,y)\
dx\ dy\right)^{-1}$ (12)

In contrast, another output of interest is the confidence collection,
explained in the subsequent subsection.

####  2.5.1 MZSM Output: Certifiable Guarantees via Confidence Collections of
Set-Based Distributions

To motivate the set-based confidence collection, notice that a single set in
the mosaic is very unlikely to hold the desired user confidence, even with a
high accuracy classifier. For example, if a classifier outputs
$\\{p(s_{j}=\mathrm{\textnormal{LOS}})\\}_{j=1}^{3}=\\{0.9,0.15,0.95\\}$ for
three satellites. Then, the leaf node matching (LOS, NLOS, LOS) would have a
probability of
$p(\lozenge^{\mathrm{\textnormal{leaf}}}_{(\mathrm{\textnormal{LOS}},\mathrm{\textnormal{NLOS}},\mathrm{\textnormal{LOS}})})=0.9\cdot
0.85\cdot 0.95=0.73$. This leaf node probability of 0.73 is much lower than a
desired confidence level of 0.95, as an example. Instead, we need to form a
collection of leaf polytopes that together achieve a user-defined confidence
level ($\gamma$, e.g., $\gamma=0.95$). We call the collection at user
confidence level $\gamma$ the confidence collection ($\mathcal{K}_{\gamma}$),
defined in Eq. (13). We use an index set
$\mathcal{L}\subseteq\\{l\\}_{l=1}^{n_{\mathrm{\textnormal{leaves}}}}$ to
select the set of leaves in the confidence collection.

$\mathcal{K}_{\gamma}=\\{\lozenge^{\mathrm{\textnormal{leaf}}}_{l}\\}_{l\in\mathcal{L}}\subseteq\mathcal{M}$
(13)

We intentionally avoid calling it a confidence interval. A confidence interval
is well defined for many unimodal distributions, such as a Gaussian
distribution. However, the symmetries of buildings in an urban plan result in
a multimodal distribution of shadow matching predictions [13, 14, 15, 39]. We
can write the process to compute the confidence collection as an optimization
problem over the index set $\mathcal{L}$ as in Eq. (14). If the user accepts
the AOI based on the self-monitor, then we can certifiably guarantee that the
user is in the confidence collection to the user’s specified confidence with
the optimization problem in Eq. (14).

$\displaystyle\min_{\mathcal{L}}$
$\displaystyle\sum_{l\in\mathcal{L}}\mathcal{O}(\lozenge^{\mathrm{\textnormal{leaf}}}_{l})$
(14) s.t.
$\displaystyle\sum_{l\in\mathcal{L}}p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l}\mid\mathcal{A})\geq\gamma$

The optimization objective $\mathcal{O}$ will depend on the user application
and the goals of the system designer. One simple choice is
$\mathcal{O}(\lozenge^{\mathrm{\textnormal{leaf}}}_{j})=1$, such that the
objective function counts the number of leaf nodes in the confidence
collection. Such an objective function seeks to find the minimal set of leaf
nodes that meets the confidence level. Intuitively, the minimal set confidence
collection obtains the level sets of the shadow matching mosaic conditioned to
the AOI. Mathematically, for the minimal set confidence collection,
$\mathcal{K}_{\gamma_{1}}\subseteq\mathcal{K}_{\gamma_{2}}$ if
$\gamma_{1}<\gamma_{2}$. In contrast, more advanced techniques can be explored
to compute confidence collections. For instance, one can use the mosaic’s tree
structure to cluster leaf polytopes or even select leaves the minimize the
area (i.e.,
$\mathcal{O}(\lozenge^{\mathrm{\textnormal{leaf}}}_{j})=\text{Area}(\lozenge^{\mathrm{\textnormal{leaf}}}_{j})$)
based on the Knapsack problem [40].

We use a greedy algorithm to find the minimal set confidence collection.
First, we take the mosaic $\mathcal{M}$ and convert it into an array. Second,
we sort the array based on descending leaf node probability
($p(\lozenge^{\mathrm{\textnormal{leaf}}}_{j}\mid\mathcal{A})$). Third, we
sequentially push leaf nodes to the confidence collection
$\mathcal{K}_{\gamma}$ from the start of the sorted array until we arrive at
the required confidence level $\gamma$. Since
$\sum_{l=1}^{n_{\mathrm{\textnormal{leaves}}}}p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l}\mid\mathcal{A})=1$,
a confidence collection $\mathcal{K}_{\gamma}$ always exists for any
confidence level $\gamma\in[0,1]$. So, this greedy algorithm is simple,
sufficient, and polynomial-time.

###  2.6 Illustrative Example of MZSM’s Memory-Efficient Tree Construction

We illustrate a simple non-physical example of MZSM in Figure 4 that includes
all overlap cases to build a clearer understanding of our proposed technique.
The AOI ($\mathcal{A}$) is a square of size $60~{}$$\mathrm{m}$$\times
60~{}$$\mathrm{m}$ oriented in cross-street and along-street directions. We
consider a total of three GNSS shadows (left side of Figure 4): (1) a
hexagonal shadow ($\mathcal{C}_{1}$) in the bottom left corner of the AOI, (2)
a hexagonal shadow in the top right corner of the AOI, and (3) a long-wide
hexagonal shadow ($\mathcal{C}_{3}$) that extends beyond the AOI. On the right
side Figure 4, we illustrate how we expand the tree as we sequentially process
each shadow. We denote branches with the no shadow overlap case with a dashed
line, the split overlap case with the solid line, and the full shadow overlap
case with a dot-dash line. The tree structure would be different with a
different satellite ordering. However, the tree always expands with one of the
three overlap cases (Proposition 1) and always yields the same four leaves
(Proposition 2). Lastly, we showcase the progression of the mosaic in the
middle of Figure 4, culminating in the tree’s output $\mathcal{M}$ (Eq. (10)).
For simplicity, we assign $p(s_{j}=\mathrm{\textnormal{NLOS}})=0.8~{}\forall
j\in\\{1,\cdots,3\\}$ and color code the leaves with the probability. Our
memory-efficient binary tree has 3 leaves at the second expansion and 4 leaves
at the third expansion. As noted in Proposition 3, these 4 non-empty leaves
fully complete the mosaic. The full binary tree would have 4 and 8 leaves at
the second and third expansions, respectively. Even in this simple example, we
already see the benefits of our MZSM’s memory-efficient binary tree
construction. From the output mosaic $\mathcal{M}$, we can formulate the PDF
(Eq. 12) or the confidence collection (Eq. 13). As noted before, a risk-aware
autonomous agent, should use the confidence collection to have certifiable
guarantees. From the confidence collection and self-monitoring probability, a
risk-aware AI decision-making agent can properly judge the uncertainty in the
scene before continuing its policy.

Figure 4: Illustrative example of a recursively-expanding binary tree with
three GNSS shadows and an LOS satellite classifier probability
$p(s_{j}=\mathrm{\textnormal{NLOS}})=0.8~{}\forall j\in\\{1,\cdots,3\\}$. The
figures on the far left column show the GNSS shadows
$\mathcal{C}_{1},\mathcal{C}_{2}$, and $\mathcal{C}_{3}$ individually overlaid
in black on the AOI $\mathcal{A}$ (gray). The middle column shows the polytope
mosaic (breakdown into different node polytopes) constructed after processing
each shadow (namely, $\mathcal{C}_{1},\mathcal{C}_{2}$, and $\mathcal{C}_{3}$
in order) based on its LOS satellite classifier probability. The node
polytopes at each tree level are colored with the probability mass function
(PMF) over the mosaic per Eq. (11). The rightmost column represents our
polytope mosaic $\mathcal{M}$ in terms of a memory-efficient binary tree
structure, where the addition of each GNSS shadow splits the node polytope
based on one of the three cases discussed in Eqs. (2)-(4), i.e., no shadow
overlap (dashed line), split overlap (solid line) and full shadow overlap
(dotted dashed line). The final layer showcases our constructed polytope
mosaic, wherein we gather the leaves from the tree and plot the leaf polytopes
that form the mosaic.

##  3 Experiments and Results

We assess the risk awareness of our proposed MZSM algorithm. We begin the
section by discussing our two baseline methods. We then detail our simulation
setup and move into the experimental results. We conduct four analyses: (1)
performance for different classifiers, (2) improvements with set-based
uncertainty models, (3) improvements with set-based confidence bounds, and (4)
self-monitoring and efficiency.

###  3.1 Baselines

This work is the first set-based risk-aware shadow matching framework, so we
cannot directly compare the set-based outputs to other methods. Specifically,
we tackle two objectives in risk-aware urban localization (Section 1): (1)
accurately modeling the uncertainty and (2) providing uncertainty bounds
across different classifier designs. GBSM is the state-of-the-art in urban
localization and serves as an intuitive baseline choice. However, we need to
isolate to what extent the risk-aware improvements from MZSM arise from
improved uncertainty modeling versus from improved uncertainty bounds.
Therefore, we use two baselines to match one baseline per objective: Gaussian
mixture models (GMMs) and set-augmented GBSM (SA-GBSM).

First, the GMM baseline isolates the first objective related to the
difficulties in handling shadow matching’s non-Gaussian uncertainties. We
provide the GMM with a precise understanding of the shadow matching
distribution and demonstrate that we cannot provide certifiable guarantees
with Gaussian error models since the Gaussian approximation does not hold.
Specifically, a GMM with one mixture component preserves the Gaussian state
error model that is often used in GBSM. For example, [27, 28, 29, 39] derive a
position and covariance estimate from the weighted average of top-scoring
cells based on an implicit single component Gaussian assumption. The multiple
mixture component cases allow us to extend beyond the current state-of-the-art
GBSM approach and demonstrate more general results about the issues of
Gaussian error models in risk-aware urban localization. We provide
implementation details of GMM in Section 3.1.1.

Second, the SA-GBSM baseline allows us to isolate the second objective related
to the uncertainty bounds across different classifiers. We base SA-GBSM on the
recent GBSM grid filter [30, 31], whereby the grid cells have different
probability mass and a selection of top-scoring grid cells define the position
update. In SA-GBSM, we augment GBSM to a set-based domain where each cell
represents a square area rather than a single point. We demonstrate how the
grid-based uncertainty bounds in SA-GBSM converge to the exact polytope set-
based uncertainty model in MZSM, though only at high grid resolutions that are
generally unpractical to compute in real-time. We provide implementation
details of SA-GBSM in Section 3.1.2.

####  3.1.1 GMM details

To provide a risk-aware analogue of GBSM that preserves the Gaussian
uncertainty model, we fit GMMs with different numbers of mixture components.
Since GMM fitting requires fitting to sampled data, we use the MZSM tree
structure to sample one million points from the exact shadow matching
distribution. We use the off-the-shelf Matlab GMM fitting function
(fitgmdist), which initializes with the $k$-means++ algorithm and then
optimizes the GMM likelihood with the iterative Expectation-Maximization
algorithm [41]. We use modest regularization to help the algorithm converge at
such large numbers of samples. We also run five replicates of Expectation-
Maximization and select the fit with the best likelihood. Generally, the
fitting process takes about one second for the 1 mixture component case, about
tens of seconds for the 2 or 3 mixture component case, about one minute for
the 4 mixture component case, and about one to four minutes for the 5 mixture
component case.

####  3.1.2 GBSM and SA-GBSM details

To provide a risk-aware analogue of GBSM that preserves the grid-based
localization, we develop SA-GBSM. SA-GBSM is nearly the same as GBSM, except
that we output a set-based localization.

We first discuss the probabilistic model shared between GBSM and SA-GBSM. GBSM
stores a building boundary look-up table for each grid cell. GBSM scores each
cell based on the probabilistic consistency between that cell and the
available satellite LOS or NLOS classification information. In particular, the
matching probability for the $j$-th satellite ($p^{\text{match}}_{j}$) is
given as Eq. (15) in works such as [23, 27]

$\displaystyle p^{\text{match}}_{j}$
$\displaystyle=p(s_{j}=\mathrm{\textnormal{LOS}}\mid
C/N_{0})p(s_{j}=\mathrm{\textnormal{LOS}}\mid
BB)+[1-p(s_{j}=\mathrm{\textnormal{LOS}}\mid
C/N_{0})][1-p(s_{j}=\mathrm{\textnormal{LOS}}\mid BB)]$ (15)

where the term $p(s_{j}=\mathrm{\textnormal{LOS}}\mid C/N_{0})$ is the
probability that the satellite is in the user’s LOS given the measurement
carrier to noise density ratio. This term can be readily switched with a
machine learning-based classifier that can use many signal features as
illustrated in [6, 20]. The second term,
$p(s_{j}=\mathrm{\textnormal{LOS}}\mid BB)$, is the probability that the
satellite is in the LOS given the building boundary look-up table. In past
GBSM works [27] where real-world experiments and dynamic obstacles were
considered, the authors set $p(s_{j}=\mathrm{\textnormal{LOS}}\mid BB)=0.85$
if the building boundary predicts that the satellite is visible and 0.2 if the
building boundary predicts that the satellite would be blocked. Intuitively,
these constants cap the range of $p^{\text{match}}_{j}$. Such caps make sense
since the GBSM building boundaries cannot be adapted in real-time and the 3D
building map does not include all the objects in the user’s environment [15].
For our current analysis, we consider a simulated set-up (as detailed in the
next section) with an ideal building model (i.e., no dynamic obstacles, etc.),
and thus set the constants as 1 and 0, respectively. The final score at each
grid is the product of the matching probabilities across satellites, namely
$\prod_{j}p^{\text{match}}_{j}$ [28, 29]. As noted in our methods (Section 2),
our MZSM scoring reflects the product of LOS/NLOS classifier probabilities
across satellites, which largely mirrors the state-of-the-art GBSM scoring
[29].

We modify the GBSM framework into SA-GBSM to treat each grid cell as a square
with side lengths matching the grid resolution following the motivation for
grid filters in recent GBSM work [30, 31]. Then, each grid cell’s square is a
set, and the squares together fully tile the AOI. Usually with GBSM, the user
must calibrate the number of top-scoring candidates based on the grid spacing,
number of satellites, and building geometry. With our set-based confidence
collections, we do not require further calibration once the designer assigns
the clustering objective (Section 2.5.1). We use the same minimal set
confidence collection objective for both MZSM and SA-GBSM (Section 2.5.1).

###  3.2 Simulated Experiment Set-up

We use a publicly-available 3D building map of San Francisco [42] and we
isolate a relevant subsection of the 3D map near the Salesforce tower
(translucent pink box at the bottom left in Figure 5), which is a notable
landmark in San Francisco. Salesforce tower sits at the intersection of
Mission and Fremont Streets, where eight large buildings (translucent pink
boxes in Figure 5) form an urban canyon. We consider the ground truth of
receiver position at $(0,-18)~{}$m in local map coordinates (black circle
marker in Figure 5) and the AOI ($\mathcal{A}$) to be a square of
$120~{}$m$\times 120~{}$m (black box outline in Figure 5) while excluding the
building footprints of six buildings that lie within the AOI. Using the same
software simulator set-up [43, 44] as explained in Figure 9 of our prior work
[32], we generate $14$ GPS satellites. Based on the relevant subsection of the
3D building map and the chosen ground truth of receiver location, 10 among a
total of 14 satellites are expected to exhibit NLOS characteristics. Next, we
utilize our prior work [32] to perform the following steps: a) convert the
entire 3D building map from their standard vertex representation to
constrained zonotope representation, b) compute the 2D GNSS shadow (a
constrained zonotope) for each satellite/building pair (from among 14
satellites and 8 buildings) by extending the 3D building constrained zonotope
in the “shadow direction” from the satellite to the building, c) convert the
shadow of each satellite/building pair from constrained zonotope to a vertex
representation that is used for representing 2D polytopes, and d) concatenate
the vertex-represented regions of 2D GNSS shadows across all the buildings for
each satellite, and then convert them to polygons of GNSS shadow (one for each
satellite). From there, we proceed with MZSM as detailed in Section 2 and as
illustrated in Section 2.6.

###  3.3 Assessing and Predicting the Shadow Matching Mosaic for different
classifiers

(a) Classifier Posterior of 50%
(Uninformative Classifier)

(b) Classifier Posterior of 75%
(Typical $C/N_{0}$ Bayesian Classifier)

(c) Classifier Posterior of 85%
(Typical SVM Classifiers)

(d) Classifier Posterior of 99.99%
(Near-Perfect Classifier)

Figure 5: PMF of our shadow matching mosaic
$p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l}\mid\mathcal{A})$ (as defined in
Eq. (11)) for four classifier posteriors of 50%, 75%, 85% and 99.99%, where
the classifier true positive and true negative rates are assumed equivalent.
The translucent magenta boxes are the building footprints, the thick black box
outline is the AOI bounding box, and the labeled black circle is the ground
truth receiver location. Our proposed MZSM validates the presence of two
distinct modes in (d) that matches the results from an ideal classifier case
discussed in our prior work [32]. In contrast, as expected, our MZSM forms a
uniform PMF over the entire mosaic using an uninformed classifier in (a),
i.e., posterior of 50%. Comparing the classifier posterior case of 75% in (b)
to that of 85% in (c), we demonstrate a noticeably more precise expected
localization with the concentration of the probability mass to fewer leaves.

With the vast space of classifiers, it is not practical to evaluate the risk-
awareness of our proposed algorithm for each possible classifier. Instead, we
take the expected value of the shadow matching mosaic ($\mathcal{M}$, defined
in Eq. (10)) via the classifier’s true positive and negative rates.
Specifically, in expectation, the classifier’s posterior over predictions
converge to its true positive rate for LOS detection and its true negative
rate for NLOS detection, by definition of the true positive and negative
rates. So, we can input the true positive and negative rates as the node
probabilities into the tree expansion of the LOS and NLOS branches, per
Section 2. We assume that the true positive rate is equal to the true negative
rate throughout the subsequent results, i.e., the true positive rate and true
negative rate are symmetric. Consistently in the results, we shorthand this
approach as “Classifier Posterior of $p\%$” to mean that we take the expected
value of the mosaic given a classifier with an accuracy of $p\%$ and with
symmetric true positive and negative rates. Note that one can extend these
results to consider asymmetric rates without loss of generality.

The classifier posterior method to evaluate the shadow matching mosaic is an
enabling technique during system design and online AI decision making. First,
the classifier posterior method allows the system designer to assess where in
the city shadow matching is most or least informative. It also allows the
designer to rapidly test different confidence collection objectives (Eq. 14)
for different AI-driven LOS classifiers without having to run each classifier
in the loop with shadow matching. Most importantly, it allows a system
designer to properly propagate the uncertainty in AI-driven LOS classifiers
into the state uncertainty to enable downstream AI decision-making. Second,
the classifier posterior method allows an AI agent to reason about how the
shadow matching uncertainty can evolve over the course of a trajectory. We
note that our tree-based structure in concert with the set-based mosaic
uniquely enables this classifier posterior method. In GBSM, there is no
embedded structure to connect grid cells based on the available satellite
information or shadow geometry. Grid cells with similar probability may or may
not hold identical LOS/NLOS classification. In MZSM, our tree provides this
necessary structure and allows the agent to forecast the mosaic, the PMF over
the mosaic, and the confidence collections. Throughout these results we take
the role of a system design for AI agents where we seek to understand the
expected shadow matching mosaic and adjust the design of an AI agent with
insights from MZSM.

We assess our proposed algorithm’s performance under four values of classifier
posteriors (Figure 5). The lower bound on classifiers is an uninformative
classier with a 50% classifier posterior. An uninformative classifier produces
a uniform distribution over the scene and does not locate the user. The upper
bound on classifiers is a near-perfect classifier with a very high classifier
posterior, such as 99.99%. The near-perfect classifier produces a more precise
user location and mirrors our past work in [32] where we used an ideal
classifier with 100% classifier posterior. In between, we have more typical
classifiers as detailed in [6, 20]. In particular, we consider a classifier
posterior of 75% that reflects a simple one-dimensional Bayesian classifier
with $C/N_{0}$, whose accuracy varies between 56% to 96% depending on the
scene and exterior building material. Lastly, we consider an 85% classifier
posterior to model the state-of-the-art support vector machine (SVM)
approaches, which use various satellite signal characteristics to consistently
achieve 80% to 90% classification accuracy. A system designer should use their
AI-driven classifier’s performance at their city of interest in more tailored
analyses.

Our analyses in Figure 5 demonstrate that the jump from 75% to 85% classifier
accuracy allows for a noticeably more precise expected localization by
concentrating the probability mass to the fewer leaves. The two disjoint modes
(most clearly seen as the red polytopes in the near-perfect classifier case)
represent the same leaf. In this scenario, there is not sufficient information
from shadow matching to provide a more precise localization than the near-
perfect classifier. Our past work [32] provides further insights into this
disjoint set behavior. With MZSM, the 75%, 85%, and 99.99% classifier
posterior cases locate the correct street and correct side of the street, in
expectation.

###  3.4 Improved risk-awareness in uncertainty models via MZSM compared to
GMMs

(a) Mosaic PDF
($f_{\mathcal{M}}$)

(b) GMM PDF ($k=1$)
($f_{\text{GMM},\,{1}}$)

(c) GMM PDF ($k=2$)
($f_{\text{GMM},\,{2}}$)

(d) GMM PDF ($k=3$)
($f_{\text{GMM},\,{3}}$)

(e) GMM PDF ($k=4$)
($f_{\text{GMM},\,{4}}$)

(f) GMM PDF ($k=5$)
($f_{\text{GMM},\,{5}}$)

Figure 6: GMM fits to the MZSM mosaic across a number of mixture components
($k=1,2,3,4,5$) for a 85% classifier posterior, which matches a typical SVM
classifier. Note that, while Figure 5(c) also discusses a classifier posterior
of 85%, it plots PMF over the mosaic
$p(\lozenge^{\mathrm{\textnormal{leaf}}}_{l}\mid\mathcal{A})$ from Eq. (11).
In contrast, we showcase the PDF over the mosaic $f_{\mathcal{M}}(x,y)$
conditioned on the AOI, per Eq. (12), in a). The thick white box outline marks
the finite support of the mosaic, which is the AOI bounding box. A GMM has
infinite support and crosses over the mosaic’s finite support, however, we
show the truncated GMM to match the AOI information in the mosaic. In (c), the
GMMs correctly locate the two modes of the shadow matching mosaic for $k\geq
2$. In comparison to (a), we demonstrate that by increasing the mixture
components from $k=1$ in (b) to $k=5$ in (f) we cannot match the preciseness
of our proposed MZSM with a GMM. Even with the 5 mixture component GMM
($f_{\text{GMM},\,{5}}$) overestimates the probability in the adjacent streets
and alleys while underestimating the mosaic at the polytope edges.

Gaussian error models are highly popular due to modeling simplicity but can be
far too restrictive of an assumption. Per figure 5, the shadow matching mosaic
is highly non-Gaussian. In Figure 6, we show the results of fitting a GMM to
the PDF of the shadow matching mosaic ($f_{\mathcal{M}}$, defined in Eq. (12))
at a classifier posterior of 85%. As illustrated in Figure 5, the shadow
matching mosaic ($\mathcal{M}$) visually appears to have two modes. The two
visual modes are more apparent as the red polytopes in the near-perfect
classifier case. With two mixture components, the optimizer fits one mixture
component of roughly equal weight (52% and 48%) around each of the two
visually identified modes. With four mixture components, it adds low-weighted
mixture components (weights of 2% and 8%) around each visually identified
mode. The mode with a 2% weight has a wide covariance to match the spread
across the entire street and the mode with an 8% weight is more centered
between the two visually identified modes to match the spread across the
intersection. With five mixture components, it separates the mode nearer to
the center of the AOI into three mixture components with weights of 22%, 22%,
and 2%. Even then, the Gaussian mixture model with 3 to 5 mixture components
appears quite similar to the two mixture component Gaussian mixture model.
None of the Gaussian mixture models match the PDF of the shadow matching
mosaic ($f_{\mathcal{M}}$). Since performance plateaus, we only entertain GMMs
with up to 5 mixture components.

We can quantify to what extent the GMM approximates the shadow matching mosaic
at a given classifier posterior as a percent error. Specifically, we can
numerically integrate the error between the $k$-mixture-component GMM
probability density ($f_{\text{GMM},\,k}$) versus the probability density
($f_{\mathcal{M}}$) of shadow matching mosaic ($\mathcal{M}$) (Equation 17).

$\displaystyle\delta_{\%}$ $\displaystyle=\left(\iint_{\text{AOI
Box}}f_{\mathcal{M}}(x,y)\ dx\ dy\right)^{-1}\iint_{\text{AOI
Box}}|f_{\text{GMM},\,k}(x,y)-f_{\mathcal{M}}(x,y)|\ dx\ dy$ (16)
$\displaystyle\delta_{\%}$ $\displaystyle=\iint_{\text{AOI
Box}}|f_{\text{GMM},\,k}(x,y)-f_{\mathcal{M}}(x,y)|\ dx\ dy$ (17)

where $\delta_{\%}$ is the integrated percent error, $f_{\text{GMM},\,k}(x,y)$
is the GMM PDF for $k$ mixture components, AOI Box is the AOI bounding box,
and ($x,y$) are the cross-street and along-street dimensions, respectively. We
condition $f_{\text{GMM},\,k}(x,y)$ on the AOI bounding box to integrate to
unity in finite integration bounds. Simply, this truncates the GMM to the AOI
bounding box as seen in Figure 6. By construction, $\delta_{\%}$ is bounded to
$0\%\leq\delta_{\%}\leq 200\%$ (equivalently $0\leq\delta_{\%}\leq 2$). When
$\delta_{\%}>100\%$, the distributions are significantly anti-aligned such
that one is high when the other is low. When $\delta_{\%}<100\%$, the
distributions are roughly aligned such that the GMM generally captures the
modes of the MZSM distribution. A $\delta_{\%}>50\%$ would be a poor
approximation, a $10\%<\delta_{\%}<50\%$ would be a decent approximation, and
$\delta_{\%}<10\%$ would be a good approximation.

######  Proposition 4.

The integrated percent error ($\delta_{\%}$) is bounded to
$0\leq\delta_{\%}\leq 2$.

###### Proof.

See D. ∎

Figure 7: Integrated Percent Error of GMM approximations of the MZSM mosaic to
evaluate the GMM approximation accuracy. We compare the GMM approximation
across mixture components ranging from 1 to 5 and across five classifier
posteriors of 50%, 75%, 85%, 95%, and 99.99%. Across all classifier
posteriors, the GMM approximation greatly improves from a 1 mixture component
model to a 2 mixture component model, but the percent error remains consistent
thereafter. None of the GMM models reach a percent error below 50% and
therefore do not represent a decent approximation of our shadow matching
mosaic.

In Figure 7, we evaluate the integrated percent error ($\delta_{\%}$) for up
to 5 mixture components and classifier posteriors of 50%, 75%, 85%, 95%, and
99.99%, where we displayed the GMM and MZSM comparison qualitatively for the
85% classifier posterior case in Figure 6. As observed qualitatively in Figure
6, the one component Gaussian is largely anti-aligned with the shadow matching
mosaic. The GMM is high at its mean is roughly between the two modes of
mosaic, but low elsewhere. Correspondingly, $\delta_{\%}>90\%$ for all the
classifiers (Figure 7). We see an improvement moving from a single Gaussian
mixture component to two Gaussian mixture components with $\delta_{\%}$
dropping from between 90% and 160% to between 60% and 110%. However,
performance plateaus thereafter. A good approximation should have
$\delta_{\%}<10\%$ or better. Even with 5 mixture components, the Gaussian
mixture model does not even hit $\delta_{\%}=50\%$. Therefore, a Gaussian
mixture model can identify all the modes and also provide a rough estimate of
the mean of the modes of the shadow matching distribution. However, a Gaussian
mixture model fails to provide a reasonable approximation of the overall
shadow matching mosaic and is thus, not suitable for risk-aware urban
localization with shadow matching. As an immediate corollary, a Gaussian
distribution or a bank of Gaussian distributions in any further filtering
framework (i.e., a Kalman filter) is not suitable for use with risk-aware
shadow matching. An AI decision-making agent using a Gaussian state model with
shadow matching is vulnerable to underestimating or overestimating the
uncertainty in the scene. A set-based uncertainty model is necessary to equip
the autonomous agent with the accurate information needed for risk-awareness.

###  3.5 Comparison analysis of Set-Based Confidence Bounds with MZSM and SA-
GBSM

In this subsection, we assess the benefits of set-based techniques to compute
confidence bounds. As detailed in Section 2.5.1, we take a collection of sets
such that the summation of probabilities associated with the sets meets the
user confidence level. We call this collection the confidence collection,
whose outer edge is the boundary of the union of all the sets in the
collection (purple outline in Figures 8 and 9).

(a) 68% Confidence Collection ($1\sigma$)
$\mathcal{K}_{\gamma={0.68}}$

(b) 90% Confidence Collection
$\mathcal{K}_{\gamma={0.90}}$

(c) 95% Confidence Collection ($2\sigma$)
$\mathcal{K}_{\gamma={0.95}}$

(d) 99.7% Confidence Collection ($3\sigma$)
$\mathcal{K}_{\gamma={0.997}}$

Figure 8: Confidence collections at different user-defined confidence levels
$\gamma=0.68,0.9,0.95,\text{ and }0.997$ for our shadow matching mosaic
($\mathcal{M}$) produced using a classifier posterior of 85%. We showcase the
selected set of leaf polytopes in the confidence collection, per Eq. (13),
where the colors of different leaf polytopes indicate the AOI-conditioned PMF
(similar to that shown in Figure 5(c)). In each case, we denote the union of
all the sets in the collection, whose outer edge is shown by a thick purple
line. We demonstrate that our proposed MZSM is confident in the side of street
up to 95% confidence ($2\sigma$) and confident in the street up to 99.7%
confidence ($3\sigma$).

In Figure 8, we illustrate four confidence collections for the shadow matching
mosaic. At a 68% ($1\sigma$) user confidence level, the confidence collection
is disjoint in the position domain, clustered around the two visually
identified modes. By the 90% user confidence level, the confidence collection
is no longer disjoint. As the user confidence level increases to 95%
($2\sigma$) and 99.7% ($3\sigma$) the confidence collection includes more leaf
polytopes and expands. For a classifier with 85% accuracy, we find that we can
determine the side of the street up to a $2\sigma$ confidence collection and
we can determine the street up to a $3\sigma$ confidence collection.

(a) Full GBSM Grid
3 $\mathrm{m}$ resolution

(b) 95% Confidence Collection $\mathcal{K}_{\gamma={0.95}}$
for GBSM at 3 $\mathrm{m}$ resolution

(c) Full GBSM Grid
5 $\mathrm{m}$ resolution

(d) 95% Confidence Collection $\mathcal{K}_{\gamma={0.95}}$
for GBSM at 5 $\mathrm{m}$ resolution

Figure 9: GBSM (as detailed in [27, 28, 29]) for a classifier posterior of
85%, which matches a typical SVM classifier. The subfigures a) and c) denote
the GBSM for a grid resolution of 3 m and 5 m, respectively, while the
subfigures b) and d) denote the selected grid points that form the confidence
collection for GBSM at a user-defined confidence level $\gamma=0.95$,
respectively. In b) and d), we denote the outer edge of all the selected grid
points by a thick purple line. As the grid resolution refines, the associated
confidence collection more closely resembles our proposed MZSM’s performance
(Figure 8(c)).

With SA-GBSM, we can also compute confidence collections, as detailed in
Section 3.1.2 The precision of the SA-GBSM confidence collection depends on
the grid resolution. In Figure 9, we illustrate the entire grid-based shadow
matching output and the 95% confidence collection for grid resolutions of 3
and 5 meters. We use 3 and 5 meter resolutions to match prior GBSM work (e.g.,
[14, 21, 28]). Qualitatively, the grid-based 95% confidence collections have a
similar shape to the MZSM 95% confidence collection in Figure 8. However, the
grid-based 95% confidence collections are three disjoint sets and do not fully
match any sharp edges in the exact MZSM 95% confidence collection.
Quantitatively, we can assess the grid-based approximation with the
“intersection over union” (IOU) metric (also known as the Jaccard index) [45].
For sets $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, the IOU metric is

$\text{IOU}(\mathcal{D}_{1},\mathcal{D}_{2})=\frac{\text{Area}(\mathcal{D}_{1}\cap\mathcal{D}_{2})}{\text{Area}(\mathcal{D}_{1}\cup\mathcal{D}_{2})}$
(18)

where the IOU metric is unity if the sets $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$ are equaland zero if the sets are completely disjoint. We
evaluate the IOU metric of the confidence collections produced with MZSM
(i.e., the exact shadow matching mosaic $\mathcal{M}$) with respect to the
confidence collections produced with SA-GBSM (i.e., the grid-based
approximation of the shadow matching mosaic) at a variety of specified user
confidence levels and classifier posteriors. Mathematically, $\mathcal{D}_{1}$
and $\mathcal{D}_{2}$ in Eq. (18) match the union of all the sets in the
confidence collection (visually, the purple outer edges in Figures 8 and 9).
We then determine to what extent the grid-based approximation converges to the
shadow matching mosaic as we refine the grid resolution.

(a) Typical $C/N_{0}$ Bayesian Classifier

(b) Typical Support Vector Machine Classifier

(c) Top Performing Classifiers

(d) Near-Perfect Classifier

Figure 10: Convergence of confidence collection (based on IOU metric in
Eq.(18)) as the grid resolution of SA-GBSM is refined from $30~{}$m to
$3~{}$m. We analyze the variation of the IOU metric as the classifier
posterior is varied from 75% to 99.99% and user-defined confidence ($\gamma$)
from 68% to 99.7%. We observe that the SA-GBSM confidence collections converge
to the MZSM confidence collections as we refine the grid resolution from 30 m
to 3 m. However, SA-GBSM more readily converges to MZSM (i.e., the slope is
shallower) at lower classifier posteriors (e.g., 75% classifier posterior in
(a)) and higher user confidence levels (e.g., $\gamma=99.7\%$).

In Figure 10, we compare convergence in grid resolution across 4 classifier
posteriors largely matching Figure 5. Note that a 50% classifier posterior is
uninformativeso, we replace it with a 95% classifier posterior, which
corresponds with a top-performing SVM classifier in a favorable environment
[6, 20]. While past GBSM works use a grid resolution of roughly 3-5 meters,
the building boundaries of a 3-5 meter resolution grid cannot be updated
sufficiently fast for real-time applications. We, therefore, illustrate
sparser grid resolutions of 10 $\mathrm{m}$ and 30 $\mathrm{m}$ that can be
updated in less time, though at degraded precision. Overall, SA-GBSM
confidence collections converge to the MZSM confidence collections as we
refine the grid resolution from 30 $\mathrm{m}$ to 3 $\mathrm{m}$ (Figure 10)
going from IOU metrics of 0.4 or less at 30 $\mathrm{m}$ to between 0.4 and
0.8 at 3 $\mathrm{m}$. The confidence collections cover a larger area at lower
classifier posteriors and higher user confidence levels, so a less refined
grid resolution better covers the larger area. Correspondingly, the IOU metric
is higher and the convergence slope is shallower at lower classifier
posteriors and higher user confidence levels. Using our novel MZSM approach,
we achieve the performance of highly refined grid resolutions in SA-GBSM while
avoiding discretization trade-offs and ensuring compute efficiency.

###  3.6 Self-Monitoring and Efficiency with Tree-based Construction

We demonstrate the benefits of the tree-based structure: how we self-monitor
the root of the tree based on the satellite consistency, how we track the tree
expansion, and how we quantify the computational load of MZSM.

Figure 11: Our MZSM’s estimate of AOI violation probability $p_{\emptyset}$,
per Eq. (9), for different classifier posteriors, when the true receiver
position is in the root AOI. Less accurate classifiers (i.e., classifier
posterior of 50%) provide a less confident assertion that the receiver is in
AOI, thus the AOI violation probability is a high value $\approx 1$. In
contrast, more accurate classifiers (i.e., classifier posterior of 99.99%)
provide a highly confident assertion, thus AOI violation probability is very
low $\approx 0$.

The AOI sets the scale for the MZSM localization task. Ideally, we want a
smaller AOI without compromising risk-aware localization. By design, the tree
estimates the probability that the GNSS shadow information violates the
assumption that the receiver is in the AOI, which is the root of the MZSM tree
(Section 2). In Figure 11, we demonstrate how this probability changes with
the classifier posterior. Specifically, we focus on the current scenario where
the true receiver location is in the AOI. A highly accurate receiver (i.e.,
classifier posterior of 99.99%) concentrates the probability on the leaves
that match the true signal designation as LOS or NLOS (such as seen in Figure
5(d)). In such a case, we accept that the receiver location is in the AOI with
the self-monitor logging the AOI violation probability at 0.001. A less
accurate receiver (i.e., classifier posterior of 50%) spreads the probability
mass over many leaves (such as seen in Figure 5(a)), and thus, the self-
monitor alerts that the AOI violation probability is 0.988. Practically, the
user should consider larger AOIs in expectation when employing a less accurate
classifier. As illustrated in Figure 11, the relation is nonlinear. At
moderate classifier accuracies, a small improvement in classifier accuracy
yields a noticeably lower AOI violation probability such that the tree
provides a noticeably more confident assertion that the receiver is in the
root AOI. Importantly, our AOI self-monitoring adapts to the classifier and
environment, thereby allowing risk awareness without being dependent on a
specific classifier.

(a) Increase in the number of leaves as the tree expands.

(b) Sorted leaf areas in the final layer of the tree.

Figure 12: Leaf analysis as the tree expands to a depth of 14 satellites. In
(a), we track the number of leaves (blue line with circular markers) at each
iteration of the tree expansion. We fit a quadratic (orange line) to
illustrate the algorithm’s quadratic complexity. In (b), we sort the leaves by
area and plot them sequentially. Only 102 leaves of 196 total have an area
larger than 5 ${\mathrm{m}}^{2}$ and only 86 leaves have an area larger than
10 ${\mathrm{m}}^{2}$ (where 10 ${\mathrm{m}}^{2}$ is just above the area of a
sedan’s footprint).

Our tree tracks its expansion with new satellite information. In Figure 12, we
compute the number of leaves after each expansion of the tree. As noted in
Section 2, the tree leaves expand quadratically with the tree depth. We
experimentally demonstrate this complexity (Figure 12(a)). The number of
leaves that the tree expansion step generates depends on how much complete
overlap there is between the new shadow and any of the past shadows. If the
new shadow is similar to any of the past shadows, many leaf nodes will be
completely in the new shadow or completely out of the new shadow. Therefore,
the tree will generate fewer new leaves (matching the no shadow overlap and
full shadow overlap cases discussed in Section 2). In contrast, a satellite
that provides unique geometry will have a shadow that intersects many nodes
and produce more new leaves (matching the split overlap case discussed in
Section 2). In our simulation, the 7th and 11th satellites are blocked
throughout most of the AOI. The shadows from these two satellites are both
similar to the 3rd satellite. Correspondingly, the 7th and 11th satellites did
not generate many new leaves at tree layers 7 and 11, respectively. In
contrast, the 9th satellite has a unique shadow and generated many more leaves
at tree layer 9. Nevertheless, the tree expansion is still quadratic. As the
number of shadows increases, some leaf nodes will have a very low area (Figure
12(b)). By the 14th layer, only a few leaves have large areas and many leaves
have small areas. For example, only 102 leaves of the 196 total leaves have an
area larger than 5 ${\mathrm{m}}^{2}$, which is roughly half the area of a
sedan footprint.

Figure 13: Timing results to complete each new layer of the tree. As the
number of layers increases, the algorithm must traverse more nodes. However,
we maintain the quadratic complexity of the tree expansion in the timing
results. Some satellites are very informative and require many intersections
yielding larger times to complete the layer (e.g., layer number 9). Other
satellites are not as informative and predominantly update the probabilities
(e.g., layer number 11).

Figure 13 illustrates the time to complete each tree expansion as the tree
depth increases, wherein our proposed algorithm produces the full tree and the
mosaic in a low median time of $0.63~{}$$\mathrm{s}$. As observed in Figure
12(a), the shadows from the 7th and 11th satellites were not as informative as
other satellites, such as the 9th satellite. We see this trend repeat in the
timing analysis. When expanding to layers 7 and 11, we identify the redundant
information with the 3rd satellite very early in the tree. So, we compute
fewer intersections and instead quickly recursively propagate the LOS and NLOS
scores to the leaves (Algorithm 1). When expanding to layer 9, the satellite
shadow is highly informative and we can split many leaf nodes. Therefore, we
produce many more leaves and we incur higher timing costs in computing the
intersections. By design, our algorithm spends more time on the more
informative signals that improve the risk-aware localization and much less
time on signals that have less impact on the risk-aware localization.

##  4 Conclusion

We proposed a novel Mosaic Zonotope Shadow Matching (MZSM) algorithm that
advanced shadow matching to enable risk-aware 3DMA-GNSS localization. We
specifically targeted two risk-aware objectives: (1) formulate precise models
of the shadow matching uncertainty in a computationally efficient manner while
tractably scaling with increasing numbers of available signals and (2) provide
guarantees on the uncertainty bounds of the urban localization solution across
a vast space of AI-driven LOS classifier designs. At a high level, we met
these objectives through set-based representations of the user location, set-
based uncertainty bounds on the user location, and efficient tree-based
construction methods for organizing the set-based information. We started with
the exact set-based representations of the GNSS shadows from our prior work
[32], where our set-based shadow extraction method avoids the trade-off
between discretization and computational load. We expanded upon this set-based
shadow matching paradigm to produce a shadow matching mosaic. We recursively
expanded a binary tree that exhaustively considered all possible LOS and NLOS
misclassifications. As the tree expands, we refined the AOI (i.e., the coarse
initial estimate) into smaller polytopes with associated probability masses.
We used special cases of overlap between the GNSS shadow and tree node
polytope to design a memory-efficient binary tree that grows quadratically
with the number of satellites rather than the exponential growth of a full
binary tree. With a probability mass function over the tree’s leaves, we
incorporated self-monitoring into the tree to check the AOI violation
probability. Finally, we derived confidence collections that provide
guarantees on the uncertainty bounds of the urban localization solution.

We validated that our MZSM algorithm meets our risk-aware objectives via our
high-fidelity experiments. First, we demonstrated that we can assess the
expectation of the shadow matching mosaic with data on the LOS/NLOS
classifier’s true positive and true negative rates (Figure 5). We discussed
how our classifier posterior method can enable a system designer to ascertain
the benefits of different classifiers and enable AI decision-making agents to
forecast how the shadow matching mosaic could evolve along a trajectory. We
employed our classifier posterior method in the role of a system designer to
perform a sensitivity analysis of the expectation of our shadow matching
mosaic across classifier accuracy. Second, we demonstrated that the Gaussian
uncertainty models, such as Gaussian mixture models (GMMs), can capture the
multimodality in the shadow matching mosaic, but they cannot provide
sufficient approximation of the uncertainty in the exact shadow matching
mosaic (Figure 6). Quantitatively, we observed that the integrated percent
error between the GMMs and the mosaic is 50% or worse (Figure 7), which is far
too high for a risk-aware approximation. Such high approximation error impedes
downstream AI agents from making risk-aware decisions and would yield
vulnerability in safety-critical applications. Hence, our exact set-based
model of the uncertainty advances shadow matching to incorporate risk
awareness. Third, we demonstrated that grid-based shadow matching approaches
can provide meaningful approximations of the uncertainty at high grid
resolutions (Figure 8 and 9). Quantitatively, we used the intersection over
union (IOU) metric to compute the similarity in the overlap between
uncertainty bounds from our MZSM method and existing work. We found that the
IOU metric is between 0.5 and 0.8 for most confidence levels for typical AI-
driven classifiers with a 3-meter grid resolution (Figure 10). Reaching better
overlap between the grid-based approximation and exact set-based mosaic would
require much finer grid resolutions. However, such fine grid resolution would
not be computationally efficient. Hence, our exact set-based model is better
suited to provide a risk-aware understanding of the uncertainty in a
computationally efficient matter than existing methods. Lastly, we
demonstrated the self-monitoring, tractability, and timing of our algorithm.
We showed how our algorithm accurately estimates the self-monitoring violation
probability to be a low value of 0.001 for a near-perfect classifier and a
high value of 0.988 for an uninformative classifier (Figure 11). Furthermore,
we validated that our tree grows quadratically with the number of satellites
with a low quadratic coefficient of 1.2 (Figure 12). We completed the tree in
a median of 0.63 seconds with the 3rd quartile at 0.65 seconds for a case with
14 shadows (Figure 13). Altogether, we succeeded in meeting our objectives and
advancing shadow matching for risk-aware 3DMA-GNSS that is well suited for the
safety of autonomous systems in harsh urban environments. At a macroscopic
level, we provide an instance where we formally connect the uncertainty in AI-
tools for localization to the state uncertainty information needed for AI
decision-making agents. At least in the domain of 3DMA-GNSS, we show the value
of certifiable and risk-aware localization in enabling risk-aware autonomous
systems.

##  5 Acknowledgments

This material is based upon work supported by the National Science Foundation
under Grant No. DGE-1656518. The authors would like to thank Adyasha Mohanty
for reviewing drafts of this paper. We would further like to thank the rest of
the Stanford Navigation and Autonomous Vehicles laboratory, Shreyas Kousik,
and Erik Herrera for insightful discussions and feedback.

##  Appendix A : Proof of Proposition 1

###### Proof.

First, we prove that the cases are distinct assuming that
$\lozenge_{j-1,n}\neq\emptyset$ and $\mathcal{C}_{j}\neq\emptyset$. The ‘No
Shadow Overlap’ case defined in Eq. (2) is the only case that satisfies
$\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\emptyset$, so it is distinct from the
rest. Similarly, the ‘Full Shadow Overlap’ case definied in Eq. (4) is the
only case that satisfies
$\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\lozenge_{j-1,n}$, so it is distinct from
the rest. Lastly, the ‘Split Overlap’ case is necessarily distinct from the
‘No Shadow Overlap’ case via
$\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\emptyset$ and necessarily distinct
from the ‘Full Shadow Overlap’ case via
$\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\lozenge_{j-1,n}$, per Eq. (3).
Therefore, the three cases are distinct, as desired. Note that if
$\lozenge_{j-1,n}=\emptyset$ or $\mathcal{C}_{j}=\emptyset$, then we have
degeneracy.

Second, we prove that the cases are complete, i.e., that the union over all
cases covers all possibilities. Without loss of generality (since unions are
commutative), consider the union of the ‘No Shadow Overlap’ case and the
‘Split Overlap’ case:

$\displaystyle\big{(}\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\emptyset\big{)}\cup\big{(}(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\emptyset)\cap(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\lozenge_{j-1,n})\big{)}$
$\displaystyle=$ (A.1)
$\displaystyle\big{(}(\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\emptyset)\cup(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\emptyset)\big{)}\cap\big{(}(\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\emptyset)\cup(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\lozenge_{j-1,n})\big{)}$
$\displaystyle=$ (A.2)
$\displaystyle\big{(}(\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\emptyset)\cup(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\lozenge_{j-1,n})\big{)}$
$\displaystyle=$ (A.3)
$\displaystyle(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\lozenge_{j-1,n})$ (A.4)

Now take the union between the above and the ‘Full Shadow Overlap’ case.
Simply,

$(\lozenge_{j-1,n}\cap\mathcal{C}_{j}\neq\lozenge_{j-1,n})\cup(\lozenge_{j-1,n}\cap\mathcal{C}_{j}=\lozenge_{j-1,n})$
(A.5)

is a tautology. Therefore, the cases are complete, as desired. ∎

##  Appendix B : Proof of Proposition 2

###### Proof.

Our binary search tree consists of branching based on two operations: set
subtraction and set intersection. Set subtraction is an instance of set
intersection. Namely, for given sets $A$ and $B$, $A-B=A\cap B^{\mathsf{c}}$,
where $B^{\mathsf{c}}$ denotes the complement of $B$. Furthermore, the
intersection of sets is commutative: for sets $A$, $B$, and $C$, $A\cap B\cap
C=A\cap C\cap B=C\cap A\cap B=\cdots$, etc. Vertex-represented polytopes are
closed under set intersection, so these intersection commutativity properties
holds for polytopes.

As seen in Eq. (B.1), we can unravel the recursion to illustrate that any node
polytope $\lozenge_{j,n}$ is computed via a chain of set intersections.

$\lozenge_{j,n}=\mathcal{A}\cap\mathcal{C}_{1}^{\
(\mathsf{c})}\cap\mathcal{C}_{2}^{\
(\mathsf{c})}\cap\cdots\cap\mathcal{C}_{j-1}^{\
(\mathsf{c})}\cap\mathcal{C}_{j}^{\ (\mathsf{c})}$ (B.1)

where (c) toggles complement for the LOS or NLOS cases per Eqs. (5) and (6).
Unfurling the recursive structure means that a position node polytope defined
via above Eq. (B.1) has a probability score given by Eq. (B.2).

$p(\lozenge_{j})=p(\mathcal{A})\cdot
p(s_{1}=(\mathrm{\textnormal{LOS}})^{(\mathsf{c})})\cdot
p(s_{2}=(\mathrm{\textnormal{LOS}})^{(\mathsf{c})})\cdots
p(s_{j-1}=(\mathrm{\textnormal{LOS}})^{(\mathsf{c})})\cdot
p(s_{j}=(\mathrm{\textnormal{LOS}})^{(\mathsf{c})})$ (B.2)

Suppose we have $k$ total GNSS shadows. By intersection commutativity, we can
equivalently construct the leaf node by any of the following intersection
chains

$\displaystyle\lozenge_{\mathrm{\textnormal{leaf}}}$
$\displaystyle=\mathcal{A}\cap\mathcal{C}_{1}^{\
(\mathsf{c})}\cap\mathcal{C}_{2}^{\
(\mathsf{c})}\cap\cdots\cap\mathcal{C}_{(k-1)}^{\
(\mathsf{c})}\cap\mathcal{C}_{k}^{\ (\mathsf{c})}$ (B.3)
$\displaystyle=\mathcal{A}\cap\mathcal{C}_{k}^{\
(\mathsf{c})}\cap\mathcal{C}_{1}^{\
(\mathsf{c})}\cap\cdots\cap\mathcal{C}_{k-1}^{\
(\mathsf{c})}\cap\mathcal{C}_{2}^{\ (\mathsf{c})}$ (B.4)
$\displaystyle=\mathcal{C}_{k}^{\ (\mathsf{c})}\cap\mathcal{C}_{1}^{\
(\mathsf{c})}\cap\cdots\cap\mathcal{C}_{k-1}^{\
(\mathsf{c})}\cap\mathcal{A}\cap\mathcal{C}_{2}^{\ (\mathsf{c})}$ (B.5)
$\displaystyle=\text{any shuffling thereof}$

This holds for all the leaves. Consider the fully expanded $2^{n}$ binary tree
at depth $n_{\mathrm{\textnormal{sat}}}$. Each leaf has a unique branching
path and, therefore, a unique identification via the complement toggle ((c))
and AOI. When we change the ordering of the satellites or AOI, we simply
change the order in which we perform the intersections. We maintain a
bijection between the original ordering and the new order such that there is a
unique mapping between a leaf in the original ordering and the equivalent leaf
in the new ordering. Therefore, the fully expanded $2^{n}$ binary tree is
commutative. As an immediate corollary, if any leaf is empty in the fully
expanded $2^{n}$ binary tree, it is empty in all satellite orderings.

Next, consider our memory-efficient binary tree. Via the same bijection as the
fully expanded $2^{n}$ binary tree, there is a unique mapping between a (non-
empty) leaf in the original ordering and the equivalent (also non-empty) leaf
in the new ordering. The only difference is that the identification is
compressed. In particular, the ‘No Shadow Overlap’ or ‘Full Shadow Overlap’
cases compress the intersection chain. In the ‘No Shadow Overlap’ case, the
$j$-th GNSS shadow ($\mathcal{C}_{j}$) does not intersect with the prior chain
of length $j$ (where we include the AOI in the chain). Hence, the LOS Child
takes the compressed identification of the parent. Combining Eq. (2), Eq. (5),
and Eq. (B.1) yields

$\displaystyle\lozenge^{\mathrm{\textnormal{LOS child}}}_{j,m}$
$\displaystyle=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}\cap\mathcal{C}_{j}^{\mathsf{c}}$
(B.6) $\displaystyle\lozenge^{\mathrm{\textnormal{LOS child}}}_{j,m}$
$\displaystyle=\big{(}\mathcal{A}\cap\mathcal{C}_{1}^{\
(\mathsf{c})}\cap\mathcal{C}_{2}^{\
(\mathsf{c})}\cap\cdots\cap\mathcal{C}_{j-1}^{\
(\mathsf{c})}\big{)}\cap\mathcal{C}_{j}^{\mathsf{c}}$ (B.7)
$\displaystyle\lozenge^{\mathrm{\textnormal{LOS child}}}_{j,m}$
$\displaystyle=\big{(}\mathcal{A}\cap\mathcal{C}_{1}^{\
(\mathsf{c})}\cap\mathcal{C}_{2}^{\
(\mathsf{c})}\cap\cdots\cap\mathcal{C}_{j-1}^{\
(\mathsf{c})}\big{)}=\lozenge^{\mathrm{\textnormal{parent}}}_{j-1,n}$ (B.8)

The equivalent holds for the ‘Full Shadow Overlap’ case and the NLOS Child.
Importantly, the compressed version results in the same polytope at the leaf.
Therefore, we have a bijection between the (non-empty) leaves in the memory-
efficient binary tree and the equivalent (non-empty) leaves in the fully
expanded $2^{n}$ binary tree. We then also have a bijection between leaves
across different satellite ordering for the memory-efficient binary tree.
Therefore, since vertex-represented polytopes are closed under set
intersection and set intersection is commutative, the leaf polytopes generated
from our binary tree are the same for any ordering of the AOI and GNSS
shadows. ∎

##  Appendix C : Proof of Proposition 3

###### Proof.

We prove the statement via induction. The base case is simply the tree’s root
(i.e., $0$-th node at $0$-th tree level), which matches the initial AOI. We
consider the following inductive hypothesis: if the union of the leaf
polytopes at the $j$-th tree expansion matches the AOI polytope $\mathcal{A}$,
then the union of the leaf polytopes at the $(j+1)$-th tree expansion also
matches the AOI polytope $\mathcal{A}$. By branching the $j$-th tree level to
the next $(j+1)$-th tree level, we split each parent node polytope into a
maximum of two children nodes: a LOS child (Eq. 5) and an NLOS child (Eq. 6).
Irrespective of the overlap case, we observe the union of the children node
polytopes equal the parent node polytope, which is shown as follows:

$\displaystyle\text{No Shadow Overlap:}\quad\lozenge^{\mathrm{\textnormal{LOS
child}}}_{j+1,m}\cup\lozenge^{\mathrm{\textnormal{NLOS child}}}_{j+1,l}$
$\displaystyle=\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}\cup\emptyset=\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}$
(C.1) $\displaystyle\text{Split
Overlap:}\quad\lozenge^{\mathrm{\textnormal{LOS
child}}}_{j+1,m}\cup\lozenge^{\mathrm{\textnormal{NLOS child}}}_{j+1,l}$
$\displaystyle=\Big{(}\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}\cap\mathcal{C}_{j+1}^{c}\Big{)}\cup\Big{(}\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}\cap\mathcal{C}_{j+1}\Big{)}=\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}$
(C.2) $\displaystyle\text{Full Shadow
Overlap:}\quad\lozenge^{\mathrm{\textnormal{LOS
child}}}_{j+1,m}\cup\lozenge^{\mathrm{\textnormal{NLOS child}}}_{j+1,l}$
$\displaystyle=\emptyset\cup\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}=\lozenge^{\mathrm{\textnormal{parent}}}_{j,n}$
(C.3)

Note that we use the distribute law to analyze the split overlap case. Namely,
for polytopes $A$, $B$, and $C$, $(A\cap B)\cup(A\cap C)=A\cap(B\cup C)$. This
result, that the union of the children polytopes matches the parent polytope,
holds for any tree expansion since these overlap cases are distinct and
complete (Propostition 1). Moreover, the empty leaves do not geometrically
contribute to completing the mosaic. By commutativity of union, if the union
of the leaf polytopes at the $j$-th tree expansion matches the AOI polytope
$\mathcal{A}$ and the union of the children polytopes matches the parent
polytopes for all tree expansions, then it must be the case that the union of
the leaf polytopes at the $(j+1)$-th tree expansion also matches the AOI
polytope $\mathcal{A}$. Therefore, by induction, the union of the leaf
polytopes matches the area of interest for any number of shadows
(equivalently, any number of satellites). ∎

##  Appendix D : Proof of Proposition 4

###### Proof.

First, we explicitly write how we condition $f_{\text{GMM},\,k}(x,y)$ on the
AOI:

$\displaystyle f_{\text{GMM},\,k}(x,y)$
$\displaystyle=\widetilde{f_{\text{GMM},\,k}}(x,y)\left(\iint_{\text{AOI}}\widetilde{f_{\text{GMM},\,k}}(x,y)\
dx\ dy\right)^{-1}$ (D.1)

where $\widetilde{f_{\text{GMM},\,k}}(x,y)$ is the Gaussian mixture model with
infinite support. We do the same with $f_{\mathcal{M}}(x,y)$ and
$\widetilde{f_{\mathcal{M}}}(x,y)$ via Eq. (12).

Now, we prove $0\leq\delta_{\%}\leq 2$. For the lower bound, notice
$f_{\text{GMM},\,k}(x,y)>0$ and $f_{\mathcal{M}}(x,y)\geq 0$. Therefore,
$\delta_{\%}\geq 0$. Equivalently, this lower bound holds from the reverse
triangle identity. For the upper bound, we use the triangle inequality
directly,

$\displaystyle|f_{\text{GMM},\,k}(x,y)-f_{\mathcal{M}}(x,y)|$
$\displaystyle\leq|f_{\text{GMM},\,k}(x,y)|+|-f_{\mathcal{M}}(x,y)|$ (D.2)
$\displaystyle\iint_{\text{AOI}}|f_{\text{GMM},\,k}(x,y)-f_{\mathcal{M}}(x,y)|\
dx\ dy$ $\displaystyle\leq\iint_{\text{AOI}}|f_{\text{GMM},\,k}(x,y)|\ dx\
dy+\iint_{\text{AOI}}|f_{\mathcal{M}}(x,y)|\ dx\ dy$ (D.3)
$\displaystyle\iint_{\text{AOI}}|f_{\text{GMM},\,k}(x,y)-f_{\mathcal{M}}(x,y)|\
dx\ dy$ $\displaystyle\leq\iint_{\text{AOI}}f_{\text{GMM},\,k}(x,y)\ dx\
dy+\iint_{\text{AOI}}f_{\mathcal{M}}(x,y)\ dx\ dy$ (D.4)
$\displaystyle\iint_{\text{AOI}}|f_{\text{GMM},\,k}(x,y)-f_{\mathcal{M}}(x,y)|\
dx\ dy$ $\displaystyle\leq 1+1$ (D.5) $\displaystyle\delta_{\%}$
$\displaystyle\leq 2$ (D.6)

Hence, $0\leq\delta_{\%}\leq 2$, as desired. ∎

## References

  * [1] T. Zhang, Q. Li, C.-s. Zhang, H.-w. Liang, P. Li, T.-m. Wang, S. Li, Y.-l. Zhu, C. Wu, Current trends in the development of intelligent unmanned autonomous systems, Frontiers of information technology & electronic engineering 18 (1) (2017) 68–85.  doi:10.1631/fitee.1601650.
  * [2] Y. Ma, Z. Wang, H. Yang, L. Yang, Artificial intelligence applications in the development of autonomous vehicles: a survey, IEEE/CAA Journal of Automatica Sinica 7 (2) (2020) 315–329.  doi:10.1109/JAS.2020.1003021.
  * [3] X. Xin, J. Jiang, Y. Zou, A review of visual-based localization, in: Proceedings of the 2019 International Conference on Robotics, Intelligent Control and Artificial Intelligence, 2019, pp. 94–105.
  * [4] G. Anderson, A. Burnheimer, V. Cicirello, D. Dorsey, C. Dugan, I. Howley, M. Kam, J. Kopena, R. Lass, K. Malfettone, et al., The AI technologies of the philadelphia area urban wireless network testbed, in: The Proceedings of the Twentieth National Conference on Artificial Intelligence, Vol. 4, 2005, pp. 1674–1675.
  * [5] R. Huang, J. Pedoeem, C. Chen, YOLO-LITE: A real-time object detection algorithm optimized for non-GPU computers, in: 2018 IEEE International Conference on Big Data (Big Data), IEEE, 2018, pp. 2503–2510.  doi:10.1109/bigdata.2018.8621865.
URL https://doi.org/10.1109%2Fbigdata.2018.8621865

  * [6] H. Xu, G. Zhang, B. Xu, L.-T. Hsu, GNSS Shadow Matching based on Intelligent LOS/NLOS Classifier, in: The 16th World Congress of the International Association of Institutes of Navigation (IAIN), Chiba, Japan, 2018, pp. 1–7.
  * [7] S. Ruan, C. Long, J. Bao, C. Li, Z. Yu, R. Li, Y. Liang, T. He, Y. Zheng, Learning to generate maps from trajectories, in: Proceedings of the AAAI conference on artificial intelligence, Vol. 34, 2020, pp. 890–897.
  * [8] M. Abdar, F. Pourpanah, S. Hussain, D. Rezazadegan, L. Liu, M. Ghavamzadeh, P. Fieguth, X. Cao, A. Khosravi, U. R. Acharya, V. Makarenkov, S. Nahavandi, A review of uncertainty quantification in deep learning: Techniques, applications and challenges, Information Fusion 76 (2021) 243–297.  doi:10.1016/j.inffus.2021.05.008.
URL https://doi.org/10.1016%2Fj.inffus.2021.05.008

  * [9] P. Misra, P. Enge, Global Positioning System: signals, measurements and performance (2nd ed. ), Ganga-Jamuna Press doi:10.1007/bf02106512.
  * [10] N. Zhu, J. Marais, D. Bétaille, M. Berbineau, GNSS position integrity in urban environments: A review of literature, IEEE Transactions on Intelligent Transportation Systems 19 (9) (2018) 2762–2778.  doi:10.1109/tits.2017.2766768.
  * [11] P. D. Groves, Shadow matching: A new GNSS positioning technique for urban canyons, NAVIGATION: Journal of the Institute of Navigation 64 (3) (2011) 417–430.  doi:10.1017/s0373463311000087.
  * [12] S. Miura, L.-T. Hsu, F. Chen, S. Kamijo, GPS error correction with pseudorange evaluation using three-dimensional maps, IEEE Transactions on Intelligent Transportation Systems 16 (6) (2015) 3104–3115.  doi:10.1109/tits.2015.2432122.
  * [13] P. D. Groves, M. Adjrad, Performance assessment of 3d-mapping–aided GNSS part 1: Algorithms, user equipment, and review, NAVIGATION: Journal of the Institute of Navigation 66 (2) (2019) 341–362.  doi:10.1002/navi.288.
  * [14] P. D. Groves, L. Wang, M. Adjrad, C. Ellul, GNSS shadow matching: The challenges ahead, in: Proceedings of the 28th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2015), 2015, pp. 2421–2443.
  * [15] M. Adjrad, P. D. Groves, J. C. Quick, C. Ellul, Performance assessment of 3D-mapping-aided GNSS part 2: Environment and mapping, NAVIGATION 66 (2) (2019) 363–383.  doi:10.1002/navi.289.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/navi.289

  * [16] L.-C. Chen, T.-A. Teo, C.-Y. Kuo, J.-Y. Rau, Shaping polyhedral buildings by the fusion of vector maps and LiDAR point clouds, Photogrammetric Engineering & Remote Sensing 74 (9) (2008) 1147–1157.  doi:10.14358/pers.73.9.1147.
  * [17] X. Chen, Y. Shao, D. He, W. Yu, Reflection objects sensing and localization with GNSS multipath signals, in: China Satellite Navigation Conference (CSNC 2021) Proceedings, Springer, 2021, pp. 203–214.  doi:10.1007/978-981-16-3146-7_20.
  * [18] C. Dostal, K. Yamafune, Photogrammetric texture mapping: A method for increasing the fidelity of 3D models of cultural heritage materials, Journal of Archaeological Science: Reports 18 (2018) 430–436.  doi:10.1016/j.jasrep.2018.01.024.
  * [19] X. Bai, W. Wen, L.-T. Hsu, Using sky-pointing fish-eye camera and LiDAR to aid GNSS single-point positioning in urban canyons, IET Intelligent Transport Systems 14 (8) (2020) 908–914.  doi:10.1049/iet-its.2019.0587.
  * [20] H. Xu, A. Angrisano, S. Gaglione, L.-T. Hsu, Machine learning based LOS/NLOS classifier and robust estimator for GNSS shadow matching, Satellite Navigation 1 (2020) 1–12.  doi:10.1186/s43020-020-00016-w.
  * [21] L. Wang, P. D. Groves, M. K. Ziebart, Urban Positioning on a Smartphone: Real-time Shadow Matching Using GNSS and 3D City Models, in: 26th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, TN, 2013, p. 14.
  * [22] L. Wang, P. D. Groves, M. K. Ziebart, GNSS shadow matching: Improving urban positioning accuracy using a 3D city model with optimized visibility scoring scheme, NAVIGATION: Journal of the Institute of Navigation 60 (3) (2013) 195–207.  doi:10.1002/navi.38.
  * [23] H.-F. Ng, G. Zhang, Y. Luo, L.-T. Hsu, Urban positioning: 3D mapping-aided GNSS using dual-frequency pseudorange measurements from smartphones, NAVIGATION 68 (4) (2021) 727–749.  doi:10.1002/navi.448.
  * [24] F. Van Diggelen, J. Wang, Improving urban GPS accuracy for your app (2020).
URL https://android-developers.googleblog.com/2020/12/improving-urban-gps-
accuracy-for-your.html

  * [25] T. Suzuki, N. Kubo, GNSS positioning with multipath simulation using 3D surface model in urban canyon, in: Proceedings of the 25th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2012), 2012, pp. 438–447.
  * [26] S. Miura, S. Hisaka, S. Kamijo, GPS multipath detection and rectification using 3D maps, in: 16th International IEEE Conference on Intelligent Transportation Systems (ITSC 2013), IEEE, 2013, pp. 1528–1534.  doi:10.1109/itsc.2013.6728447.
  * [27] L. Wang, P. D. Groves, M. K. Ziebart, Smartphone shadow matching for better cross-street GNSS positioning in urban environments, NAVIGATION: Journal of the Institute of Navigation 68 (3) (2015) 411–433.  doi:10.1017/s0373463314000836.
  * [28] M. Adjrad, P. Groves, Intelligent urban positioning using shadow matching and GNSS ranging aided by 3D mapping, in: Proceedings of the 29th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2016), Portland, OR, 2016, pp. 534–553.  doi:10.33012/2016.14845.
  * [29] M. Adjrad, P. D. Groves, Intelligent urban positioning: Integration of shadow matching with 3D-mapping-aided GNSS ranging, The Journal of Navigation 71 (1) (2018) 1–20.  doi:10.1017/s0373463317000509.
  * [30] P. D. Groves, M. Adjrad, Multi-epoch 3D mapping aided GNSS using a grid filter, in: 31st International Technical Meeting of The Satellite Division of the Institute of Navigation (ION GNSS+ 2018), Miami, Florida, 2018, pp. 3335–3356.  doi:10.33012/2018.16070.
  * [31] Q. Zhong, P. D. Groves, Multi-epoch 3D-mapping-aided positioning using bayesian filtering techniques, in: Proceedings of the 34th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2021), 2021, pp. 195–225.  doi:10.33012/2021.17894.
  * [32] S. Bhamidipati, S. Kousik, G. Gao, Set-valued shadow matching using zonotopes, in: Proceedings of the 34th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2021), 2021, pp. 2373–2390.  doi:10.33012/2021.17933.
  * [33] D. Shi, T. Chen, L. Shi, On set-valued Kalman filtering and its application to event-based state estimation, IEEE Transactions on Automatic Control 60 (5) (2014) 1275–1290.  doi:10.1109/tac.2014.2370472.
  * [34] J. K. Scott, D. M. Raimondo, G. R. Marseglia, R. D. Braatz, Constrained zonotopes: A new tool for set-based estimation and fault detection, Automatica 69 (2016) 126–136.  doi:10.1016/j.automatica.2016.02.036.
  * [35] M. Althoff, J. M. Dolan, Online verification of automated road vehicles using reachability analysis, IEEE Transactions on Robotics 30 (4) (2014) 903–918, doi:10.1109/tro.2014.2312453.  doi:10.1109/tro.2014.2312453.
  * [36] P. L. Combettes, The foundations of set theoretic estimation, Proceedings of the IEEE 81 (2) (1993) 182–208.  doi:10.1109/5.214546.
  * [37] M. Adjrad, P. D. Groves, Enhancing least squares GNSS positioning with 3D mapping without accurate prior knowledge, NAVIGATION, Journal of the Institute of Navigation 64 (1) (2017) 75–91.  doi:10.1002/navi.178.
  * [38] OEIS Foundation Inc., Central polygonal numbers (the Lazy Caterer’s sequence), Entry A000124, The On-Line Encyclopedia of Integer Sequences.
URL https://oeis.org/A000124

  * [39] L. Wang, Kinematic GNSS Shadow Matching Using a Particle Filter, in: Proceedings of the 27th International Technical Meeting of the ION Satellite Division, Tampa, FL, 2014, p. 13.
  * [40] S. Martello, P. Toth, Knapsack Problems: Algorithms and Computer Implementations, Wiley-Interscience Series in Discrete Mathematics and Optimization, J. Wiley & Sons, Chichester ; New York, 1990.  doi:10.1016/0377-2217(91)90159-s.
  * [41] MathWorks, Fit Gaussian mixture model to data - MATLAB, [Online] Available: https://www.mathworks.com/help/stats/fitgmdist.html (2022).
  * [42] R. György, 3DWarehouse San Francisco, [Online] Available: https://3dwarehouse.sketchup.com/model/4ad4796d-8102-4bdd-9bfc-5442dee9facf/San-Francisco?hl=en&login=true (2018 (Accessed as of Sept, 2021)).
  * [43] K. Paul, Soft GNSS, [Online] Available: https://github.com/kristianpaul/SoftGNSS ((Accessed as of June, 2020)).
  * [44] T. Ebinuma, GPS-SDR-SIM, [Online] Available: https://github.com/osqzss/gps-sdr-sim ((Accessed as of June, 2020)).
  * [45] H. Rezatofighi, N. Tsoi, J. Gwak, A. Sadeghian, I. Reid, S. Savarese, Generalized intersection over union: A metric and a loss for bounding box regression, in: 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), IEEE, Long Beach, CA, USA, 2019, pp. 658–666.  doi:10.1109/CVPR.2019.00075.