# Spectrum Surveying: Active Radio Map
Estimation with Autonomous UAVs

Raju Shrestha1, Daniel Romero1, and Sundeep Prabhakar Chepuri2
1Department of Information and Communication Technology,
University of Agder, Norway.
2Department of Electrical Communication Engineering,
Indian Institute of Science, India Research funded by the Research Council of
Norway (IKTPLUSS grant 280835) and Department of Science and Technology,
India. Emails: {raju.shrestha<EMAIL_ADDRESS><EMAIL_ADDRESS>Parts of this work were presented at the IEEE International Workshop on
Machine Learning for Signal Processing 2020 [1].

###### Abstract

Radio maps find numerous applications in wireless communications and mobile
robotics tasks, including resource allocation, interference coordination, and
mission planning. Although numerous techniques have been proposed to construct
radio maps from spatially distributed measurements, the locations of such
measurements are assumed predetermined beforehand. In contrast, this paper
proposes _spectrum surveying_ , where a mobile robot such as an unmanned
aerial vehicle (UAV) collects measurements at a set of locations that are
actively selected to obtain high-quality map estimates in a short surveying
time. This is performed in two steps. First, two novel algorithms, a model-
based online Bayesian estimator and a data-driven deep learning algorithm, are
devised for updating a map estimate and an uncertainty metric that indicates
the informativeness of measurements at each possible location. These
algorithms offer complementary benefits and feature constant complexity per
measurement. Second, the uncertainty metric is used to plan the trajectory of
the UAV to gather measurements at the most informative locations. To overcome
the combinatorial complexity of this problem, a dynamic programming approach
is proposed to obtain lists of waypoints through areas of large uncertainty in
linear time. Numerical experiments conducted on a realistic dataset confirm
that the proposed scheme constructs accurate radio maps quickly.

###### Index Terms:

Radio maps, spectrum cartography, UAV communications, deep learning,
trajectory planning.

##  I Introduction

Radio maps are functions that provide a certain channel metric, such as
received signal strength, power spectral density, or channel gain, across a
given geographical area. Radio maps are commonly used in tasks such as network
planning, interference coordination, power control, spectrum management,
resource allocation, hand-off procedure design, dynamic spectrum access, and
cognitive radio [2, 3], to name a few. Besides, radio maps are gaining
popularity for autonomous unmanned aerial vehicle (UAV) communications in
tasks such as mission planning [4] or optimal relay deployment in UAV-assisted
networks; see, e.g. [5, 6, 7, 8]. All these applications require methods for
accurate radio map construction.

A large number of approaches have emerged to address the radio map estimation
problem. They rely on measurements acquired by spatially distributed sensors,
possibly integrated into user equipment such as mobile phones, to construct
radio maps by means of some form of spatial interpolation. Schemes to
construct _power_ maps, which provide the received signal strength across
space, have been developed using kriging [9, 10, 11, 1], dictionary learning
[12], sparse Bayesian learning [13, 14, 15], and matrix completion [16]. Power
spectral density (PSD) maps can be estimated using kernel-based learning [17,
18, 19], sparse learning [18], and tensor completion [20, 21]. The problem of
estimating channel-gain maps has been addressed in [22, 23, 24]. All the
aforementioned approaches are based on interpolation algorithms or rely on
modeling propagation phenomena [25]. However, the actual radio propagation
environment is often complex and inappropriate modeling will generally lead to
poor estimation performance. To bypass this problem, [26, 27, 28, 29] proposed
data-driven approaches where a deep neural network (DNN) is used to learn the
underlying propagation phenomena such as shadowing, reflection, and
diffraction from a dataset of measurements.

All the preceding methods assume that the measurement locations are given and,
as a result, cannot decide where to measure next, which is necessary in
certain contexts. Specifically, certain applications such as network planning
or those involving UAV communications demand approaches to construct radio
maps by _surveying_ the channel at the area of interest with a sensor on board
a mobile robot. To this end, this work puts forth an active sensing method
where an autonomous vehicle, such as a UAV, equipped with a sensing and
communication module collects measurements at a judiciously selected set of
locations to efficiently construct a high-quality radio map. In particular,
the set of most informative measurement locations is selected to approximately
minimize the operation time required to attain a certain estimation accuracy.
Towards this purpose, suitable metrics that provide the uncertainty associated
with every spatial location are developed. Using the collected measurements,
the UAV computes this metric on-the-fly and plans a trajectory accordingly.

The contributions111The conference version of this work [1] presents the idea
of the online Bayesian algorithm for model-based uncertainty mapping and
spectrum surveying with a UAV. Relative to [1], the present paper additionally
proposes a data-driven uncertainty mapping scheme using a DNN, a new trade-off
scheme for trajectory planning, and extensive empirical validation and
comparison with existing benchmarks through a realistic dataset. of this work
are as follows:

  * C1)

An online Bayesian learning scheme is proposed for joint radio map estimation
and uncertainty mapping with constant complexity per measurement. This is
required since UAVs process measurements as they become available and update
their trajectory accordingly. The proposed algorithm constitutes an online
alternative to the popular kriging estimator [30], which is grounded on the
well-known Gudmundson shadowing model [31].

  * C2)

Complementing the benefits of the aforementioned online Bayesian estimator, a
_data-driven_ approach for jointly estimating a power map and uncertainty
metric is developed based on a DNN. As existing DNN-based radio map
estimators, the proposed algorithm learns propagation phenomena from a dataset
and, therefore, it yields a superior performance relative to the online
Bayesian estimator at the expense of increasing computational complexity.

  * C3)

A trajectory planning strategy is proposed for measurement acquisition at the
most informative locations as indicated by the uncertainty metrics from C1) or
C2). To reduce the computational complexity of such a task, a simple waypoint-
search approach based on the Bellman-Ford shortest path algorithm [32] is
presented.

The code necessary to reproduce all experiments is available at
https://github.com/uiano/spectrum_surveying_with_UAVs and a video illustrating
a spectrum surveying operation is available at https://youtu.be/r9zDr4O0Fp8.

The main novelty of this work is twofold: i) this paper proposes data-driven
uncertainty mapping in realistic propagation environments with possibly non-
Gaussian channel distribution, and ii) this paper proposes uncertainty-aware
trajectory planning for spectrum surveying with autonomous UAVs.

The rest of the paper is structured as follows: Sec. II formulates the problem
and outlines the proposed approach. Secs. III and IV respectively present
model-based and data-driven algorithms for power map estimation and
uncertainty mapping. These algorithms lay the grounds for the uncertainty-
aware trajectory planning scheme in Sec. V. The proposed scheme is empirically
validated in Sec. VI by means of data obtained from a ray-tracing simulator.
Finally, Sec. VII draws connections with related work and Sec. VIII concludes
the paper.

Notation:$|\mathcal{X}|$ denotes the cardinality of set $\mathcal{X}$. The
function $\lfloor x\rfloor$ denotes the largest integer less than or equal to
$x$. Bold uppercase letters denote matrices or tensors, bold lowercase letters
represent column vectors, and non-bold lowercase letters denote scalars.
$[\bm{x}]_{i}$ is the $i$-th entry of vector $\bm{x}$, $\text{vec}(\bm{X})$
denotes vectorization of matrix $\bm{X}$, $[\bm{X}]_{i,j}$ is the
$({{i},{j}})$-th entry of matrix $\bm{X}$, and $[\bm{Y}]_{i,j,k}$ is the
$(i,j,k)$-th entry of tensor $\bm{Y}$. The Hadamard product is represented by
$\odot$ and $\|\bm{X}\|_{F}$ refers to the Frobenius norm of $\bm{X}$.

##  II Spectrum Surveying

This section formulates the problem and outlines the proposed approach.

###  II-A Problem Formulation

Let ${\mathcal{X}}\subset\mathbb{R}^{{{{m}}}}$ be the geographical area of
interest, where ${{{m}}}$ is typically $2$ or $3$, and consider $S$
transmitters, which may correspond to cellular base stations. A power map is a
function ${r}({\bm{x}})$ which provides the value of the received power, also
known as signal strength, at every ${\bm{x}}\in{\mathcal{X}}$. A UAV equipped
with a communication module capable of measuring the received power and a
positioning system such as GPS adaptively chooses future measurement locations
based on previous measurements. The power measurements and their locations are
respectively denoted as $\\{\tilde{{r}}_{{{\tau}}}\\}_{{{\tau}}=0}^{{{t}}}$
and $\\{{\bm{x}}_{{{\tau}}}\\}_{{{\tau}}=0}^{t}$, where
$\tilde{{r}}_{{{\tau}}}={r}({\bm{x}}_{{{\tau}}})+{{z}}_{{{\tau}}}$ and
${{z}}_{{{\tau}}}$ represents measurement noise. For convenience, the
measurements up to ${{t}}+1$ will be arranged as
${{\tilde{\bm{{r}}}}}_{{{t}}}\triangleq[\tilde{{r}}_{0},\ldots,\tilde{{r}}_{t}]\in\mathbb{R}^{t+1}$
and
${\bm{X}}_{{{t}}}\triangleq[{\bm{x}}_{0},\ldots,{\bm{x}}_{t}]\in\mathbb{R}^{{{{m}}}\times(t+1)}$.
The goal is to determine ${\bm{X}}_{{{t}}}$ in order to infer an accurate
radio map as quickly as possible using the measurements collected at these
locations.

###  II-B Proposed Approach

To judiciously select the future most informative locations, one needs to
solve two sub-problems: P1) Given
$\\{({\bm{x}}_{{{\tau}}},\tilde{{r}}_{{{\tau}}}),{{\tau}}=0,...,{{t}}\\}$, the
task is to find an estimate $\hat{{r}}$ of ${r}$ and an uncertainty metric
that represents the uncertainty associated with each location. Function ${r}$
is typically referred to as the true radio map, whereas $\hat{{r}}$ is the map
estimate. An algorithm that produces $\hat{{r}}$ is termed a map estimator.
P2) Using the uncertainty metric from P1, the second problem involves planning
a trajectory for measurement collection to attain the desired estimation
accuracy as fast as possible.

The following sections will further detail the formulations of these two
problems. Secs. III and IV provide two algorithms with complementary benefits
to address P1, whereas Sec. V is concerned with P2.

##  III Model-Based Uncertainty Mapping

This section builds upon a widely-used probabilistic shadowing model to
develop an _online_ Bayesian algorithm for power map estimation and
uncertainty mapping. To find a suitable uncertainty metric, the idea is to
obtain the posterior distribution of ${r}$ given the measurements. While the
mean of such a distribution provides the _minimum mean square error_ (MMSE)
estimate of ${r}$, its variance can be used as uncertainty metric. This
contrasts with most algorithms in the literature, which just provide estimates
of ${r}$ without any associated uncertainty metric.

###  III-A Radio Map Model

To simplify the notation, consider a single transmitter with transmit power
${P}_{\text{Tx}}$. Then, the power received at ${\bm{x}}\in{\mathcal{X}}$ in
logarithmic units is given by

$\displaystyle{r}({\bm{x}})={P}_{\text{Tx}}+{{\bar{l}}}({\bm{x}})-{{\bar{{s}}}}({\bm{x}})+{{w}}({\bm{x}}),$
(1)

where ${{\bar{l}}}({\bm{x}})$ quantifies free-space path loss and antenna
gain, ${{\bar{{s}}}}({\bm{x}})$ is the loss due to shadowing, and
${{w}}({\bm{x}})$ is a gain due to small-scale fading and unmodeled effects.
The log-normal shadowing component
${{\bar{{s}}}}({\bm{x}})\sim\mathcal{N}({{{\mu}}}_{{\bar{{s}}}},{{\sigma^{2}_{{{{s}}}}}})$
follows the Gudmundson correlation model [31], which prescribes that
$\text{Cov}({{\bar{{s}}}}({\bm{x}}),{{\bar{{s}}}}({\bm{x}}^{\prime}))={{\sigma^{2}_{{{{s}}}}}}2^{-\|{\bm{x}}-{\bm{x}}^{\prime}\|/{\delta}_{{{s}}}}$,
where ${{\sigma^{2}_{{{{s}}}}}}$ is a constant and ${\delta}_{{{s}}}$ is the
distance at which the correlation decays to $1/2$. Furthermore, following
[13], ${{w}}({\bm{x}})$ is modeled as $\mathcal{N}(0,{{\sigma^{2}_{{{w}}}}})$
and assumed independent of ${{w}}({\bm{x}}^{\prime})$ and
${{\bar{{s}}}}({\bm{x}}^{\prime\prime})$ $\forall{\bm{x}}^{\prime}$,
${\bm{x}}^{\prime\prime}\in{\mathcal{X}}$ with
${\bm{x}}\neq{\bm{x}}^{\prime}$. For brevity, rewrite (1) as

$\displaystyle{r}({\bm{x}})={{l}}({\bm{x}})-{{s}}({\bm{x}})+{{w}}({\bm{x}}),$
(2)

where
${{l}}({\bm{x}})\triangleq{P}_{\text{Tx}}+{{\bar{l}}}({\bm{x}})-{{{\mu}}}_{{\bar{{s}}}}$
and
${{s}}({\bm{x}})\triangleq{{\bar{{s}}}}({\bm{x}})-{{{\mu}}}_{{\bar{{s}}}}$.
Here ${{l}}({\bm{x}})$ is assumed to be a known deterministic component since
i) ${P}_{\text{Tx}}$ and the source location can be assumed known as base
stations in contemporary cellular networks share this information with the
users and ii) ${{{\mu}}}_{{\bar{{s}}}}$ can be estimted readily from a set of
measurements.

To avoid unbounded complexity for finding estimates
$\forall{\bm{x}}\in{\mathcal{X}}$ [33, Sec. 6.4], the map and uncertainty
metrics will be evaluated at a finite set of arbitrary grid points
${\mathcal{G}}\triangleq\\{{{\bm{x}}}^{{\mathcal{G}}}_{0},\ldots,{{\bm{x}}}^{{\mathcal{G}}}_{{{{G}}}-1}\\}\subset{\mathcal{X}}$.
Thus, using (2), one can write

$\displaystyle{{\bm{{r}}}}^{{\mathcal{G}}}\triangleq[{r}({{\bm{x}}}^{{\mathcal{G}}}_{0}),\ldots,{r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{G}}}-1})]^{\top}={\bm{{{l}}}}^{{\mathcal{G}}}-{\bm{{{s}}}}^{{\mathcal{G}}}+{\bm{{{w}}}}^{{\mathcal{G}}},$
(3)

where
${\bm{{{l}}}}^{{\mathcal{G}}}\triangleq[{{l}}({{\bm{x}}}^{{\mathcal{G}}}_{0}),\ldots,{{l}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{G}}}-1})]^{\top}$,
${\bm{{{s}}}}^{{\mathcal{G}}}\triangleq[{{s}}({{\bm{x}}}^{{\mathcal{G}}}_{0}),\ldots,$
${{s}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{G}}}-1})]^{\top}$, and
${\bm{{{w}}}}^{{\mathcal{G}}}\triangleq[{{w}}({{\bm{x}}}^{{\mathcal{G}}}_{0}),\ldots,{{w}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{G}}}-1})]^{\top}$.

The UAV collects a power measurement
$\tilde{{r}}_{{{\tau}}}={r}({\bm{x}}_{{{\tau}}})+{{z}}_{{{\tau}}}$ when it is
at position ${\bm{x}}_{{{\tau}}}\in{\mathcal{X}}$, $\tau=0,1,\ldots$, where
${{z}}_{{{\tau}}}\sim\mathcal{N}(0,{{\sigma^{2}_{{{z}}}}})$ models the
measurement error and is assumed independent across ${{\tau}}$ and independent
of ${{w}}({\bm{x}})$ and ${{s}}({\bm{x}}^{\prime})$ $\forall{\bm{x}}$,
${\bm{x}}^{\prime}\in{\mathcal{X}}$. For notational purposes, let

$\displaystyle{{\tilde{\bm{{r}}}}}_{{{t}}}={\bm{{{l}}}}_{{{t}}}-{\bm{{{s}}}}_{{{t}}}+{\bm{{{w}}}}_{{{t}}}+\bm{{{z}}}_{{{t}}},$
(4)

where
${\bm{{{l}}}}_{{{t}}}\triangleq[{{l}}({\bm{x}}_{0}),\ldots,{{l}}({\bm{x}}_{{{t}}})]^{\top}$,
${\bm{{{s}}}}_{{{t}}}\triangleq[{{s}}({\bm{x}}_{0}),\ldots,{{s}}({\bm{x}}_{{{t}}})]^{\top}$,
${\bm{{{w}}}}_{{{t}}}\triangleq[{{w}}({\bm{x}}_{0}),\ldots,{{w}}({\bm{x}}_{{{t}}})]^{\top}$,
and $\bm{{{z}}}_{{{t}}}\triangleq[{{z}}_{0},\ldots,{{z}}_{{{t}}}]^{\top}$.

###  III-B Batch Bayesian Estimation for Uncertainty Mapping

To facilitate understanding, the batch version of the problem is described
before presenting the proposed online algorithm. The batch problem, commonly
referred to as _kriging_ [34], is to obtain the posterior distribution
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}},{\bm{X}}_{{{t}}})$
given ${{\tilde{\bm{{r}}}}}_{{{t}}}$ and ${\bm{X}}_{{{t}}}$. Applying a well-
known result in estimation theory [35, Ch. 10], the mean of such a posterior
distribution provides the MMSE estimate of ${r}({\bm{x}})$ at the grid points.
On the other hand, the covariance of this posterior captures the uncertainty
about the true ${{\bm{{r}}}}^{{\mathcal{G}}}$ after observing
${{\tilde{\bm{{r}}}}}_{{{t}}}$ and ${\bm{X}}_{{{t}}}$.

It is straightforward to establish that ${{\bm{{r}}}}^{{\mathcal{G}}}$ is
conditionally independent of ${{\tilde{\bm{{r}}}}}_{{{t}}}$ given
${\bm{{{s}}}}^{{\mathcal{G}}}$ using the model represented by (3) and (4). As
a result, it follows that

$\displaystyle p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})$
$\displaystyle=\int
p({{\bm{{r}}}}^{{\mathcal{G}}},{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})d{\bm{{{s}}}}^{{\mathcal{G}}}$
$\displaystyle=\int
p({{\bm{{r}}}}^{{\mathcal{G}}}|{\bm{{{s}}}}^{{\mathcal{G}}})p({\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})d{\bm{{{s}}}}^{{\mathcal{G}}},$
(5)

where ${\bm{X}}_{{{t}}}$ has been omitted for brevity. From (3) and the fact
that ${\bm{{{l}}}}^{{\mathcal{G}}}$ is deterministic, it clearly follows that
the first factor in the integrand is
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{\bm{{{s}}}}^{{\mathcal{G}}})=\mathcal{N}({{\bm{{r}}}}^{{\mathcal{G}}}|{\bm{{{l}}}}^{{\mathcal{G}}}-{\bm{{{s}}}}^{{\mathcal{G}}},{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{G}}})$.
To obtain the second factor
$p({\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})$, observe that
${\bm{{{s}}}}^{{\mathcal{G}}}$ and ${{\tilde{\bm{{r}}}}}_{{{t}}}$ are jointly
Gaussian. In particular, one can obtain the parameters of their joint
distribution $p({\bm{{{s}}}}^{{\mathcal{G}}},{{\tilde{\bm{{r}}}}}_{{{t}}})$ as
follows. Firstly, the mean vectors are
$\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}]=\bm{0}$ and
$\mathop{\mathbb{E}}[{{\tilde{\bm{{r}}}}}_{{{t}}}]={\bm{{{l}}}}_{{{t}}}$.
Next, to compute the covariance, let
$\mathop{\textrm{Cov}}[{\bm{{{s}}}}^{{\mathcal{G}}}]\triangleq{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}$
and write

$\displaystyle\mathop{\textrm{Cov}}[{\bm{{{s}}}}^{{\mathcal{G}}},{{\tilde{\bm{{r}}}}}_{{{t}}}]$
$\displaystyle=\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})^{\top}]$
$\displaystyle=\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}(-{\bm{{{s}}}}_{{{t}}}+{\bm{{{w}}}}_{{{t}}}+\bm{{{z}}}_{{{t}}})^{\top}]$
$\displaystyle=-\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}{\bm{{{s}}}}_{{{t}}}^{\top}]\triangleq-{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}$
(6)

as well as

$\displaystyle\mathop{\textrm{Cov}}[{{\tilde{\bm{{r}}}}}_{{{t}}}]$
$\displaystyle=\mathop{\mathbb{E}}[({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})^{\top}]$
$\displaystyle=\mathop{\mathbb{E}}[(-{\bm{{{s}}}}_{{{t}}}+{\bm{{{w}}}}_{{{t}}}+\bm{{{z}}}_{{{t}}})(-{\bm{{{s}}}}_{{{t}}}+{\bm{{{w}}}}_{{{t}}}+\bm{{{z}}}_{{{t}}})^{\top}]$
$\displaystyle=\mathop{\textrm{Cov}}[{\bm{{{s}}}}_{{{t}}}]+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{t}}+1}+{{\sigma^{2}_{{{z}}}}}\bm{I}_{{{t}}+1}$
$\displaystyle\triangleq{{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}}}+{{\sigma}^{2}_{{{w}},{{z}}}}\bm{I}_{{{t}}+1},$
(7)

where
${{\sigma}^{2}_{{{w}},{{z}}}}={{\sigma^{2}_{{{w}}}}}\bm{+}{{\sigma^{2}_{{{z}}}}}\bm{$}.Here,thematrices$C_
s^G,C_ s^G,  s_t$and$C_
s_t$canbeobtainedfromthecovariancefunctionintroducedinSec.~{}\ref{subsec:radiomodel}.Applying\cite[cite]{[\@@bibref{}{kay1}{}{},
Th. 10.2]}tothisjointdistribution,itfollowsthat$p( s^G—~r_t) = N( s^G—μ_
s^G—~r_t, C_ s^G—~r_t)$,with\begin{aligned}
{{\bm{\mu}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}&=\mathop{\textrm{Cov}}[{\bm{{{s}}}}^{{\mathcal{G}}},{{\tilde{\bm{{r}}}}}_{{{t}}}]{\mathop{\textrm{Cov}}}^{-1}[{{\tilde{\bm{{r}}}}}_{{{t}}}]({{\tilde{\bm{{r}}}}}_{{{t}}}-\mathop{\mathbb{E}}[{{\tilde{\bm{{r}}}}}_{{{t}}}])\\\
&=-{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}({{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}}}+{{\sigma}^{2}_{{{w}},{{z}}}}\bm{I}_{{{t}}+1})^{-1}({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})\\\
{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}&=\mathop{\textrm{Cov}}[{\bm{{{s}}}}^{{\mathcal{G}}}]-\mathop{\textrm{Cov}}[{\bm{{{s}}}}^{{\mathcal{G}}},{{\tilde{\bm{{r}}}}}_{{{t}}}]{\mathop{\textrm{Cov}}}^{-1}[{{\tilde{\bm{{r}}}}}_{{{t}}}]\mathop{\textrm{Cov}}[{{\tilde{\bm{{r}}}}}_{{{t}}},{\bm{{{s}}}}^{{\mathcal{G}}}]\\\
&={{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}-{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}({{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}}}+{{\sigma}^{2}_{{{w}},{{z}}}}\bm{I}_{{{t}}+1})^{-1}{{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}},{\bm{{{s}}}}^{{\mathcal{G}}}},\end{aligned}where$C_
s_t, s^G≜C_ s^G,
s_t^⊤$.Finally,applying\cite[cite]{[\@@bibref{}{bishop2006}{}{}, eq.
(2.115)]}toobtaintheconditionalmarginalin\eqref{eq:marginalbatch}yields$p(r^G—~r_t)
= N(r^G— μ_r^G—~r_t ,C_r^G—~r_t )$with$μ_r^G—~r_t≜l_t-μ_
s^G—~r_t$and$C_r^G—~r_t≜σ^2_wI_G+C_ s^G—~r_t
$,therebysolvingthebatchproblem.Themapestimateis,therefore,$^r^G=
μ_r^G—~r_t$and$C_r^G—~r_t$capturestheuncertaintyinthisestimate.\par\par\par$

###  III-C Online Bayesian Estimation for Uncertainty Mapping

Although the solution described in Sec. III-B could in principle be utilized
to determine the trajectory of the UAV, it suffers from a limitation: since
(III-B) involves inverting a $({{t}}+1)\times({{t}}+1)$ matrix, the complexity
per uncertainty computation will grow as more measurements become available,
eventually becoming unaffordable. Thus, it is more convenient to adopt an
_online_ approach where each new measurement is utilized to refine the
previous posterior. Specifically, consider the problem of finding
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}},{\bm{X}}_{{{t}}})$
given the previous posterior
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1},{\bm{X}}_{{{t}}-1})$
and the most recent measurement $(\tilde{{r}}_{{{t}}},{\bm{x}}_{{{t}}})$ with
a computational complexity that does not grow with ${{t}}$.

To address this problem, it is convenient to decompose
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})$ into the
previous posterior
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1})$ and a term
that depends on the current measurement but not on the previous ones. However,
it can be easily seen that such a factorization is not possible due to the
posterior correlation among measurements. To bypass this difficulty, the
central idea in the proposed online learning scheme is to use ${\mathcal{G}}$
to summarize the information of all past measurements. Mathematically, this
can be formulated as the assumption that $\tilde{{r}}_{{{t}}}$ and
${{\tilde{\bm{{r}}}}}_{{{t}}-1}$ are conditionally independent given
${{\bm{{r}}}}^{{\mathcal{G}}}$. That is, when ${{\bm{{r}}}}^{{\mathcal{G}}}$
is known, the past measurements ${{\tilde{\bm{{r}}}}}_{{{t}}-1}$ do not
provide extra information about $\tilde{{r}}_{{{t}}}$. The error that this
approximation introduces, which can be reduced by adopting a denser grid, pays
off since it enables online estimation. Besides, Sec. VI demonstrates that
this error is negligible.

From this assumption and Bayes’ rule, it follows that

$\displaystyle p({{\bm{{r}}}}^{{\mathcal{G}}}$
$\displaystyle|{{\tilde{\bm{{r}}}}}_{{{t}}})=p({{\bm{{r}}}}^{{\mathcal{G}}}|\tilde{{r}}_{{{t}}},{{\tilde{\bm{{r}}}}}_{{{t}}-1})\propto
p(\tilde{{r}}_{{{t}}},{{\tilde{\bm{{r}}}}}_{{{t}}-1}|{{\bm{{r}}}}^{{\mathcal{G}}})p({{\bm{{r}}}}^{{\mathcal{G}}})$
$\displaystyle=p(\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}})p({{\tilde{\bm{{r}}}}}_{{{t}}-1}|{{\bm{{r}}}}^{{\mathcal{G}}})p({{\bm{{r}}}}^{{\mathcal{G}}})=p({{\tilde{\bm{{r}}}}}_{{{t}}-1},{{\bm{{r}}}}^{{\mathcal{G}}})p(\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}})$
$\displaystyle=p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1})p({{\tilde{\bm{{r}}}}}_{{{t}}-1})p(\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}})\propto
p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1})p(\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}}),$
(8)

where $\propto$ denotes equality up to a scaling constant. In this context, a
constant is understood as a positive term that does not depend on
${{\bm{{r}}}}^{{\mathcal{G}}}$. The term
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1})=\mathcal{N}({{\bm{{r}}}}^{{\mathcal{G}}}|{{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}},{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}})$
is obtained from the previous iteration. To find
$p(\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}})$, note that
$\tilde{{r}}_{{{t}}}$ and ${{\bm{{r}}}}^{{\mathcal{G}}}$ are jointly Gaussian.
It follows from [35, Th. 10.2] that
$p(\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}})$ is Gaussian distributed
with parameters

$\displaystyle\mathop{\mathbb{E}}$
$\displaystyle[\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}}]=\mathop{\mathbb{E}}[\tilde{{r}}_{{{t}}}]+\mathop{\textrm{Cov}}[\tilde{{r}}_{{{t}}},{{\bm{{r}}}}^{{\mathcal{G}}}]{\mathop{\textrm{Cov}}}^{-1}[{{\bm{{r}}}}^{{\mathcal{G}}}]({{\bm{{r}}}}^{{\mathcal{G}}}-\mathop{\mathbb{E}}[{{\bm{{r}}}}^{{\mathcal{G}}}])$
$\displaystyle={{l}}({\bm{x}}_{{{t}}})+\mathop{\mathbb{E}}[(-{{s}}({\bm{x}}_{{{t}}})+{{w}}({\bm{x}}_{{{t}}})+{{z}}_{{{t}}})(-{\bm{{{s}}}}^{{\mathcal{G}}}+{\bm{{{w}}}}^{{\mathcal{G}}})^{\top}]$
$\displaystyle\quad\times{\mathop{\mathbb{E}}}^{-1}[(-{\bm{{{s}}}}^{{\mathcal{G}}}+{\bm{{{w}}}}^{{\mathcal{G}}})(-{\bm{{{s}}}}^{{\mathcal{G}}}+{\bm{{{w}}}}^{{\mathcal{G}}})^{\top}]({{\bm{{r}}}}^{{\mathcal{G}}}-{\bm{{{l}}}}^{{\mathcal{G}}})$
$\displaystyle={{\bm{a}}}_{{{t}}}^{\top}{{\bm{{r}}}}^{{\mathcal{G}}}+{{{b}}}_{{{t}}}$
(9) $\displaystyle\mathop{\textrm{Var}}$
$\displaystyle[\tilde{{r}}_{{{t}}}|{{\bm{{r}}}}^{{\mathcal{G}}}]=\mathop{\textrm{Var}}[\tilde{{r}}_{{{t}}}]-\mathop{\textrm{Cov}}[\tilde{{r}}_{{{t}}},{{\bm{{r}}}}^{{\mathcal{G}}}]{\mathop{\textrm{Cov}}}^{-1}[{{\bm{{r}}}}^{{\mathcal{G}}}]\mathop{\textrm{Cov}}[{{\bm{{r}}}}^{{\mathcal{G}}},\tilde{{r}}_{{{t}}}]$
$\displaystyle={{\sigma^{2}_{{{{s}}}}}}+{{\sigma}^{2}_{{{w}},{{z}}}}-({{\bm{C}}}_{{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\bm{C}}}_{{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}})({{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{G}}})^{-1}$
$\displaystyle\quad\times({{\bm{C}}}_{{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\bm{C}}}_{{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}})^{\top}\triangleq{{{\lambda}}}_{{{t}}},$
(10)

where the quantities
${{\bm{a}}}_{{{t}}}\triangleq({{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{G}}})^{-1}({{\bm{C}}}_{{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\bm{C}}}_{{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}})^{\top}$
and
${{{b}}}_{{{t}}}\triangleq{{l}}({\bm{x}}_{{{t}}})-({{\bm{C}}}_{{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\bm{C}}}_{{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}})({{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{G}}})^{-1}{\bm{{{l}}}}^{{\mathcal{G}}}$
have been defined along with
${{\bm{C}}}_{{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}}\triangleq\mathop{\textrm{Cov}}[{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}]$
and
${{\bm{C}}}_{{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}}\triangleq\mathop{\textrm{Cov}}[{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}]$.
Clearly, the latter contains a single non-zero entry if
${\bm{x}}_{{{t}}}\in{\mathcal{G}}$ and vanishes otherwise.

Finally, (III-C) and [33, eq. (2.116)] imply that
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})=\mathcal{N}({{\bm{{r}}}}^{{\mathcal{G}}}|{{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}},{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}})$
with

$\displaystyle{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}$
$\displaystyle=({{\bm{C}}}^{-1}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}+(1/{{{\lambda}}}_{{{t}}}){{\bm{a}}}_{{{t}}}{{\bm{a}}}_{{{t}}}^{\top})^{-1}$
$\displaystyle={{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}-\frac{{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}{{\bm{a}}}_{{{t}}}{{\bm{a}}}_{{{t}}}^{\top}{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}}{{{{\lambda}}}_{{{t}}}+{{\bm{a}}}_{{{t}}}^{\top}{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}{{\bm{a}}}_{{{t}}}}$
(11)

and

$\displaystyle{{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}$
$\displaystyle={{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}\left[\frac{\tilde{{r}}({\bm{x}}_{{{t}}})-{{{b}}}_{{{t}}}}{{{{\lambda}}}_{{{t}}}}{{\bm{a}}}_{{{t}}}+{{\bm{C}}}^{-1}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}{{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}\right].$
(12)

The sought algorithm applies the update equations (III-C) and (12) every time
a new measurement is acquired, where the initializations are given by
${{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{-1}}\triangleq{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{G}}}$
and
${{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{-1}}\triangleq{\bm{{{l}}}}^{{\mathcal{G}}}$.
Unlike the batch estimator from Sec. III-B, Equation (12) provides an
estimator
${\hat{\bm{{r}}}}^{{\mathcal{G}}}_{{t}}\triangleq{{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}$
that changes with ${{t}}$, gradually accommodating new measurement information
as it is acquired. Observe that the computational complexity is now constant
per ${{t}}$, and therefore this algorithm constitutes a practical estimator
for spectrum surveying.

##  IV Data-Driven Uncertainty Mapping

The preceding section described techniques to estimate radio maps based on a
shadowing model. Although the adopted formulation led to a practical online
estimator, the quality of the map estimate may be limited if the actual
propagation conditions do not adhere to this model. For this reason, this
section adopts a data-driven perspective where the propagation phenomena are
implicitly learnt from a training dataset of historical measurements collected
in different environments. Specifically, a deep learning approach is pursued
here to jointly obtain a radio spectral map estimate and uncertainty metric.

###  IV-A Input and Output Representation

The proposed estimator comprises a deep neural network that takes information
about the measurements at its input and returns a map estimate and an
uncertainty map at the output. This section describes how these inputs and
outputs are represented along the lines of [26, 20] to accommodate the
variable number of measurements. To this end, it is necessary to arrange the
grid points as an ${{N_{y}}\times{N_{x}}}$ rectangular grid. By doing so, each
grid point ${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ in the notation of the
previous sections can be identified with a row and a column in the rectangular
grid. Thus, the grid can be written as
${\mathcal{G}}=\\{{{\bm{x}}}^{{\mathcal{G}}}_{{{i},{j}}},{i}=0,1,\ldots,{N_{y}}-1,{j}=0,1,\ldots,{N_{x}}-1\\}$,
where
${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}={{\bm{x}}}^{{\mathcal{G}}}_{{{i},{j}}}$
if ${{{g}}}={i}{N_{y}}+{j}$. In this way, one can arrange the true received
power values at the grid points in matrix
${\bm{R}}\in\mathbb{R}^{{{N_{y}}\times{N_{x}}}}$, whose $({{i},{j}})$-th entry
is given by
$[{\bm{R}}]_{{{i},{j}}}={r}({{\bm{x}}}^{{\mathcal{G}}}_{{{i},{j}}})$.

Similarly, the measurements collected until a given time instant can be
represented by a matrix of the same dimension ${{N_{y}}\times{N_{x}}}$ by
associating each measurement location with the nearest grid point. In
particular, let
${\mathcal{T}}_{{{i},{j}}}\triangleq\\{{{\tau}}:\|{{\bm{x}}}^{{\mathcal{G}}}_{{{i},{j}}}-{\bm{x}}_{{{\tau}}}\|_{2}\leq\|{{\bm{x}}}^{{\mathcal{G}}}_{{i}^{\prime},{j}^{\prime}}-{\bm{x}}_{{{\tau}}}\|_{2}$,
$\forall{i}^{\prime},{j}^{\prime}\\}$ be a set containing the indices of the
measurement locations assigned to ${{\bm{x}}}^{{\mathcal{G}}}_{{{i},{j}}}$.
The targeted matrix representation involves associating with
${{\bm{x}}}^{{\mathcal{G}}}_{{{i},{j}}}$ the average
$\tilde{{r}}_{{{i},{j}}}^{({{t}})}=(1/|{\mathcal{T}}_{{{i},{j}}}|)\sum_{{{\tau}}\in{\mathcal{T}}_{{{i},{j}}}\cap\\{0,\ldots,{{t}}\\}}\tilde{{r}}_{{{\tau}}}$
if $|{\mathcal{T}}_{{{i},{j}}}\cap\\{0,\ldots,{{t}}\\}|>0$ and
$\tilde{{r}}_{{{i},{j}}}^{({{t}})}=0$ otherwise. The result of aggregating all
measurements in ${{\tilde{\bm{{r}}}}}_{{{t}}}$ can be expressed in matrix form
as ${\tilde{\bm{{\bm{R}}}}}_{{{t}}}\in{\mathbb{R}}^{{{N_{y}}\times{N_{x}}}}$,
where
$[{\tilde{\bm{{\bm{R}}}}}_{{{t}}}]_{{{i},{j}}}=\tilde{{r}}_{{{i},{j}}}^{({{t}})}$.
It is also convenient to introduce the notation
${\Omega}_{{{t}}}=\\{({{i},{j}}):|{\mathcal{T}}_{{{i},{j}}}\cap\\{0,\ldots,{{t}}\\}|>0\\}\subset\\{0,...,{N_{y}}-1\\}\times\\{0,...,{N_{x}}-1\\}$
to denote a set containing the indices of the grid points associated with at
least one measurement.

To accommodate general scenarios, define
${\mathcal{B}}\subset\\{0,...,{N_{y}}-1\\}\times\\{0,...,{N_{x}}-1\\}$ to be
the set of indices of the grid points that lie inside buildings or no-fly
zones. Furthermore, to encode the information in ${\Omega}_{{{t}}}$ and
${\mathcal{B}}$, a mask
${\bm{M}}_{{{t}}}\in\\{1,0,-1\\}^{{{N_{y}}\times{N_{x}}}}$ is constructed
[26], where $[{\bm{M}}_{{{t}}}]_{{{i},{j}}}=1$ if
$({{i},{j}})\in{\Omega}_{{{t}}}$, $[{\bm{M}}_{{{t}}}]_{{{i},{j}}}=-1$ if
$({{i},{j}})\in{\mathcal{B}}$, and $[{\bm{M}}_{{{t}}}]_{{{i},{j}}}=0$
otherwise. This mask is concatenated to ${\tilde{\bm{{\bm{R}}}}}_{{{t}}}$
along the third dimension to form
${\check{{\bm{R}}}}_{{{t}}}\in\mathbb{R}^{{{N_{y}}\times{N_{x}}}\times 2}$ as
the tensor input to the network.

Regarding the output, functions
${{f}^{{\bm{R}}}}:\mathbb{R}^{{{N_{y}}\times{N_{x}}}\times
2}\rightarrow\mathbb{R}^{{{N_{y}}\times{N_{x}}}}$ and
${{f}^{{\bm{V}}}}:\mathbb{R}^{{{N_{y}}\times{N_{x}}}\times
2}\rightarrow\mathbb{R}^{{{N_{y}}\times{N_{x}}}}$ will be used to denote the
radio map estimate and uncertainty metric obtained by the neural network,
which are thus encoded as ${{{N_{y}}\times{N_{x}}}}$ matrices. The map
estimate provided by the network given the measurements in
${{\tilde{\bm{{r}}}}}_{{{t}}}$ is therefore given by
${\hat{\bm{{r}}}}^{{\mathcal{G}}}_{{{t}}}={{f}^{{\bm{R}}}}({\check{{\bm{R}}}}_{{{t}}})$.
Thus, every time a new measurement is collected, it is necessary to update
${\check{{\bm{R}}}}_{{{t}}}$ and evaluate the output of the network (forward
pass).

###  IV-B Training Loss Design

The goal now is to design a suitable training loss to learn the functions
${{f}^{{\bm{R}}}}$ and ${{f}^{{\bm{V}}}}$ introduced in Sec. IV-A. Let
$\\{({\check{{\bm{R}}}}^{({d})},{\bm{R}}^{({d})})\\}_{{d}=0}^{{D}-1}$
represent pairs of measurements and true power maps collected in ${D}$
environments, where ${\check{{\bm{R}}}}^{({d})}$ and ${\bm{R}}^{({d})}$ are
respectively of the form of ${\check{{\bm{R}}}}_{{{t}}}$ and ${\bm{R}}$
described in Sec. IV-A. Subscript ${{t}}$ is omitted to simplify the notation.
Given these training pairs, one could think of adopting a negative log-
likelihood loss function to learn ${{f}^{{\bm{R}}}}$ and ${{f}^{{\bm{V}}}}$ as
the posterior mean and posterior standard deviation of ${r}$ whenever the
distribution of ${r}$ evaluated at an arbitrary set of locations is jointly
Gaussian [36]. However, this is not generally the case in practice. To see
this, Fig. 1 depicts a histogram of ${r}({\bm{x}})$ for the Rossyln dataset
[26] described in Sec. VI. It is observed that the distribution is markedly
non-Gaussian. For this reason, an alternative loss function is investigated
next.

Figure 1: Histogram of true maps from the Rossyln dataset.

The idea is to rely on the well-known relation between the MMSE estimator and
the conventional square loss. This is explained in detail next when designing
the learning approach for ${{f}^{{\bm{R}}}}$ with the purpose of laying the
grounds for ${{f}^{{\bm{V}}}}$. Specifically, recall that the MMSE estimator
of ${\bm{R}}$, i.e., the minimizer of the mean square error (MSE) between a
target ${\bm{R}}$ and its estimate, is the conditional mean
$\mathop{\mathbb{E}}[{\bm{R}}|{\check{{\bm{R}}}}]$ [35, Ch. 10]

$\displaystyle\arg\min_{{{f}^{{\bm{R}}}}}\mathop{\mathbb{E}}[\|{\bm{R}}-{{f}^{{\bm{R}}}}({\check{{\bm{R}}}})\|^{2}_{F}]=\mathop{\mathbb{E}}[{\bm{R}}|{\check{{\bm{R}}}}],$
(13)

where the right-hand side is seen as a function of ${\check{{\bm{R}}}}$ and,
to simplify the exposition, it is assumed that there are no buildings in the
map. The case with buildings is more involved since measurements cannot be
collected indoors. For this reason, it is assumed that ${r}$ is not defined
inside buildings. The case with buildings is studied later.

The left-hand side of (13) does not directly yield an estimator for two
reasons. On the one hand, the expectation in (13) is over random pairs
$({\check{{\bm{R}}}},{\bm{R}})$, but their distribution is not known. Thus,
one can replace the expectation with a sample average, which results in

$\displaystyle\arg\min_{{{f}^{{\bm{R}}}}}\mathop{\mathbb{E}}[\|{\bm{R}}-{{f}^{{\bm{R}}}}({\check{{\bm{R}}}})\|^{2}_{F}]$
$\displaystyle\approx\arg\min_{{{f}^{{\bm{R}}}}}\frac{1}{{D}}\sum_{{d}=0}^{{D}-1}\|{\bm{R}}^{({d})}-{{f}^{{\bm{R}}}}({\check{{\bm{R}}}}^{({d})})\|^{2}_{F}.$
(14)

On the other hand, the optimization in the left-hand side of (13) is over an
arbitrary function ${{f}^{{\bm{R}}}}$. Unfortunately, solving such a problem
is only tractable in very specific circumstances. Instead, one generally needs
to constrain the feasible set of functions e.g., to those parameterized by a
neural network. In this case, one can optimize over the set
$\\{{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}:{\bm{{\theta}}}\in\mathbb{R}^{P}\\}$
and find

$\displaystyle{\hat{{\bm{{\theta}}}}}$
$\displaystyle=\arg\min_{{\bm{{\theta}}}}\frac{1}{{D}}\sum_{{d}=0}^{{D}-1}\|{\bm{R}}^{({d})}-{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})})\|^{2}_{F}.$
(15)

Clearly, if the set
$\\{{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}:{\bm{{\theta}}}\in\mathbb{R}^{P}\\}$ is
sufficiently large and ${D}$ is sufficiently high, then
${{f}^{{\bm{R}}}_{{\hat{{\bm{{\theta}}}}}}}$ will be close to the MMSE
estimator.

As an uncertainty metric, one can think of estimating the posterior variance,
along the lines of Sec. III. To this end, note that the matrix collecting the
posterior variances for all grid points is given by the conditional mean of
the residual
$({\bm{R}}-\mathop{\mathbb{E}}[{\bm{R}}|{\check{{\bm{R}}}}])^{2}$, where
$(.)^{2}$ denotes entry-wise square. The key is to generalize and apply (13)
in the reverse direction as before, which yields

$\displaystyle\text{Var}[{\bm{R}}|{\check{{\bm{R}}}}]$
$\displaystyle=\mathop{\mathbb{E}}\Big{[}({\bm{R}}-\mathop{\mathbb{E}}[{\bm{R}}|{\check{{\bm{R}}}}])^{2}|{\check{{\bm{R}}}}\Big{]}$
$\displaystyle=\arg\min_{({{f}^{{\bm{V}}}})^{2}}\mathop{\mathbb{E}}\|({\bm{R}}-\mathop{\mathbb{E}}[{\bm{R}}|{\check{{\bm{R}}}}])^{2}-({{f}^{{\bm{V}}}}({\check{{\bm{R}}}}))^{2}\|^{2}_{F},$
(16)

where ${{f}^{{\bm{V}}}}$ is chosen to denote posterior standard deviation
rather than posterior variance for reasons that will become clear later. If
$\mathop{\mathbb{E}}[{\bm{R}}^{({d})}|{\check{{\bm{R}}}}^{({d})}]$ were known,
one could adopt the approach in (15) to solve

$\displaystyle{\hat{{\bm{{\theta}}}}}=\arg\min_{{\bm{{\theta}}}}\frac{1}{{D}}\sum_{{d}=0}^{{D}-1}$
$\displaystyle\|({\bm{R}}^{({d})}-\mathop{\mathbb{E}}[{\bm{R}}^{({d})}|{\check{{\bm{R}}}}^{({d})}])^{2}$
$\displaystyle-({{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})}))^{2}\|^{2}_{F}.$
(17)

In practice, though,
$\mathop{\mathbb{E}}[{\bm{R}}^{({d})}|{\check{{\bm{R}}}}^{({d})}]$ is unknown.
However, as argued earlier, it is close to
${{f}^{{\bm{R}}}_{{\hat{{\bm{{\theta}}}}}}}$, where ${\hat{{\bm{{\theta}}}}}$
is given by (15). Thus, it makes sense to replace
$\mathop{\mathbb{E}}[{\bm{R}}^{({d})}|{\check{{\bm{R}}}}^{({d})}]$ in (IV-B)
with ${{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})})$ and add
the objective in (15) as an additional term that promotes that
${{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}$ is close to
$\mathop{\mathbb{E}}[{\bm{R}}^{({d})}|{\check{{\bm{R}}}}^{({d})}]$.

This gives rise to

$\displaystyle{\hat{{\bm{{\theta}}}}}=$
$\displaystyle\arg\min_{{\bm{{\theta}}}}~{}\Big{\\{}(1-{\alpha})\frac{1}{{D}}\sum\limits_{{d}=0}^{{D}-1}\|{\bm{R}}^{({d})}-{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})})\|^{2}_{F}$
$\displaystyle+{\alpha}\frac{1}{{D}}\sum\limits_{{d}=0}^{{D}-1}\|({\bm{R}}^{({d})}-{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})}))^{2}-({{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})}))^{2}\|^{2}_{F}\Big{\\}},$
(18)

where ${\alpha}\in[0,1]$ can be adjusted to balance the trade-off between
learning ${{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}$ and
${{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}$. Observe that the second term in (IV-B)
involves squares of squares of its arguments. This may render its scale very
different from the one of the first term, which may produce numerical problems
and hinder the selection of ${\alpha}$. To alleviate such issues,
${{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}$ is trained to learn the absolute residual
$|{\bm{R}}_{{d}}-{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})})|$
rather than the squared residual
$({\bm{R}}^{({d})}-{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})}))^{2}$,
where $|\cdot|$ denotes entry-wise absolute value. This suggests training the
DNN by solving

$\displaystyle{\hat{{\bm{{\theta}}}}}=\arg\min_{{\bm{{\theta}}}}\frac{1}{{D}}\sum\limits_{{d}=0}^{{D}-1}\Big{\\{}(1-{\alpha})\|{{\Delta}_{{d},{\bm{{\theta}}}}^{{\bm{R}}}}\|^{2}_{F}+{\alpha}\|{{\Delta}_{{d},{\bm{{\theta}}}}^{{\bm{V}}}}\|^{2}_{F}\Big{\\}},$
(19)

where
${{\Delta}_{{d},{\bm{{\theta}}}}^{{\bm{R}}}}={\bm{R}}^{({d})}-{{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})})$,
and
${{\Delta}_{{d},{\bm{{\theta}}}}^{{\bm{V}}}}=|{{\Delta}_{{d},{\bm{{\theta}}}}^{{\bm{R}}}}|-{{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}^{({d})})$.

Several improvements of (19) are discussed next. Let us begin by noting that a
good map estimate must be close to the measurements at the corresponding
measurement locations, and that a reasonable uncertainty metric is expected to
yield low values near measurement locations. To promote this behavior, one may
assign a heavier weight to the residuals at the measurement locations. Let set
${\Omega}^{({d})}$ contain the indices of the measurement locations of the
${d}$-th training example and let ${\lambda}\in[0.5,1]$ represent a scaling
factor. In this case, one may consider solving

$\displaystyle{\hat{{\bm{{\theta}}}}}=\arg\min_{{\bm{{\theta}}}}\frac{1}{{D}}\sum\limits_{{d}=0}^{{D}-1}\Big{\\{}$
$\displaystyle(1-{\alpha})\|{{\Delta}_{{d},{\bm{{\theta}}}}^{{\bm{R}}}}\odot{\bm{K}}_{{d}}\|^{2}_{F}$
$\displaystyle+{\alpha}\|{{\Delta}_{{d},{\bm{{\theta}}}}^{{\bm{V}}}}\odot{\bm{K}}_{{d}}\|^{2}_{F}\Big{\\}},$
(20)

where ${\bm{K}}_{{d}}\in\mathbb{R}^{{{N_{y}}\times{N_{x}}}}$ has entries
$[{\bm{K}}_{{d}}]_{{{i},{j}}}={\lambda}$ if $({{i},{j}})\in{\Omega}^{({d})}$
and $[{\bm{K}}_{{d}}]_{{{i},{j}}}=1-{\lambda}$ otherwise. If ${\lambda}=0.5$,
the minimizer of (IV-B) boils down to the minimizer of (19). Setting
${\lambda}>0.5$ increases the focus on observed locations.

For scenarios with buildings and no-fly zones, let ${\mathcal{B}}_{{d}}$
represent the set of indices of the grid points inside buildings in the
${d}$-th training example. Since no measurements can be acquired indoors, i.e.
${\mathcal{B}}_{{d}}\cap{\Omega}^{({d})}=\emptyset$, it would make no sense to
fit the DNN to any value there. In other words, one can assign a zero weight
to the residuals inside buildings by setting
$[{\bm{K}}_{{d}}]_{{{i},{j}}}={\lambda}$ if $({{i},{j}})\in{\Omega}^{({d})}$,
$[{\bm{K}}_{{d}}]_{{{i},{j}}}=0$ if $({{i},{j}})\in{\mathcal{B}}_{{d}}$, and
$[{\bm{K}}_{{d}}]_{{{i},{j}}}=1-{\lambda}$ otherwise.

###  IV-C Network Architecture

The previous section assumed a generic family of functions
$\\{({{f}^{{\bm{R}}}_{{\bm{{\theta}}}}},{{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}),{\bm{{\theta}}}\in\mathbb{R}^{P}\\}$.
The form of these functions is next specified by developing a network
architecture suitable for uncertainty-aware radio map estimation.

At the high level, the two outputs of the DNN will be implemented here using a
separate subnetwork with dedicated parameters. Specifically, let these
subnetworks be denoted by ${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}$ and
${{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}$, whose weights are collected
respectively in ${\bm{{\theta}}}_{{\bm{R}}}$ and ${\bm{{\theta}}}_{{\bm{V}}}$.
Subnetwork ${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}$ takes
${\check{{\bm{R}}}}$ as its input and returns an estimate of ${\bm{R}}$ as its
output. The adopted architecture for this subnetwork is a convolutional
autoencoder given its well-documented merits in radio map estimation [26].

An autoencoder is a type of neural network
${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}$ that can be expressed as [37,
Ch. 14]:
${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}({\check{{\bm{R}}}})={\bar{D}}_{{{\bm{{\theta}}}_{{\bar{D}}}}}({\bar{E}}_{{{\bm{{\theta}}}_{{\bar{E}}}}}({\check{{\bm{R}}}}))$,
where functions ${\bar{E}}_{{{\bm{{\theta}}}_{{\bar{E}}}}}$ and
${\bar{D}}_{{{\bm{{\theta}}}_{{\bar{D}}}}}$ are known respectively as encoder
and decoder and ${{\bm{{\theta}}}_{{\bar{E}}}}$ and
${{\bm{{\theta}}}_{{\bar{D}}}}$ denote their associated weight vectors. The
output of the encoder
${\bm{\zeta}}={\bar{E}}_{{{\bm{{\theta}}}_{{\bar{E}}}}}({\check{{\bm{R}}}})\in\mathbb{R}^{N_{{\bm{\zeta}}}}$,
which is known as the code or vector of latent variables, is of a (possibly
much) smaller dimension $N_{{\bm{\zeta}}}$ than the dimension of
${\check{{\bm{R}}}}$. Training the autoencoder forces the encoder to condense
or summarize the information of ${\check{{\bm{R}}}}$ into just
$N_{{\bm{\zeta}}}$ variables. This is useful to exploit the prior information
that the input ${\check{{\bm{R}}}}$ lies on a low-dimensional manifold
embedded in a high dimensional space. This is the case of radio maps, as
established in [26].

The reason for implementing an autoencoder using convolutional layers [38] is
the ability of the latter to efficiently learn and exploit spatial
correlation, as reflected in their so-called _translational equivariance_
property. As detailed in [26], this is highly beneficial in the case of radio
map estimation, where the signal strength is substantially correlated at
nearby locations.

Having presented the implementation of ${{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}$,
the next step is to implement ${{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}$. To this
end, it is useful to think of which information must be fed at the input of
${{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}$. On the one hand, the
information that ${\check{{\bm{R}}}}$ conveys about the measurement locations
is intuitively relevant to determine the uncertainty metric since points that
lie far away from the measurement locations are expected to be assigned a
greater uncertainty. On the other hand, the map estimate obtained by
${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}$ may also be useful for
uncertainty metric computation since it provides information about the spatial
variability of the map. For this reason, the subnetwork
${{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}$ makes use of both
${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}({\check{{\bm{R}}}})$ and
${\check{{\bm{R}}}}$, which leads to a cascaded architecture where
${{f}^{{\bm{R}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}})={{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}({\check{{\bm{R}}}})$
and
${{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}})={{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}({{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}({\check{{\bm{R}}}}),{\check{{\bm{R}}}})$.

Fig. 2 depicts the high-level architecture of the proposed deep convolutional
neural network, which will be termed deep radio map and uncertainty estimator
(DRUE) hereafter. Subnetwork ${{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}$
is also implemented as a convolutional autoencoder since it is empirically
observed that this yields a satisfactory performance; see Sec. VI.

Figure 2: A high level architecture of a deep radio map and uncertainty
estimator (DRUE). Figure 3: A detailed view of DRUE based on an autoencoder
architecture. DRUE has two subnetworks
${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}$ and
${{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}$, where each subnetwork is a
fully convolutional autoencoder with an identical structure.

Implementation details of DRUE are shown in Fig. 3. The encoder comprises
mostly 2D convolutional and max pooling layers whereas the decoder comprises
2D convolutional transpose layers, also known as a deconvolutional layers, and
up-sampling layers. Leaky rectified linear units (Leaky ReLUs) [39] are
utilized as activations in all convolutional and deconvolutional layers except
in the output layer of ${{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}$, which
utilizes an exponential [40] activation to guarantee the non-negativity of the
uncertainty. Each convolutional and deconvolutional layer has a kernel size of
$4\times 4$, a stride of $1$, and same padding. Max pooling layers have a pool
size of 2 and a stride of 2, whereas upsampling layers have an upsampling
factor of 2.

##  V Uncertainty-Aware Trajectory Planning

The uncertainty metrics provided by the estimators in Secs. III and IV are
used in this section to develop a path planning scheme for measurement
collection with a UAV.

The uncertainty in the target map at grid point
${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ after observing
${{\tilde{\bm{{r}}}}}_{{{t}}}$ will be denoted by
${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})\in\mathbb{R}_{+}$. In the case
of the Bayesian estimators from Sec. III,
${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ will be given by

$\displaystyle{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})=\big{[}{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}\big{]}_{{{{g}}},{{{g}}}}$
(21)

if there is a single transmitter. When there are multiple transmitters, the
values of the right-hand side of (21) for all transmitters can be averaged to
obtain a single value of ${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ per
${{{g}}}$. For the estimator in Sec. IV,
${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ is given by

$\displaystyle{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})=\big{[}\text{vec}({{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}_{{{t}}}))\big{]}_{{{{g}}}},$
(22)

where ${\check{{\bm{R}}}}_{{{t}}}$ is obtained from
${{\tilde{\bm{{r}}}}}_{{{t}}}$ as described in Sec. IV-A. In this case, the
estimator directly provides the aggregate uncertainty for all transmitters
and, therefore, no averaging is required.

Finally, as a performance metric, it is convenient to define the total
uncertainty of the map after observing ${{\tilde{\bm{{r}}}}}_{{{t}}}$ as the
spatial average of the point-wise uncertainty values. Since the uncertainty
metric is only meaningful outside buildings or no-fly zones (cf. Sec. IV-B),
this average becomes

$\displaystyle{{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})\triangleq\frac{1}{|{\bar{\mathcal{B}}}|}\sum_{{{{g}}}\in{\bar{\mathcal{B}}}}{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}),$
(23)

where ${\bar{\mathcal{B}}}$ is the complement of the set ${\mathcal{B}}$ that
contains the indices of the grid points inside buildings or no-fly zones; see
Sec. IV-A.

###  V-A Trajectory Planning

Ideally, one would wish to plan a trajectory for measurement collection such
that the error in the map estimate decreases as quickly as possible. However,
knowing the error requires knowing the true map, which is not possible in
practice. Instead, one may use the uncertainty metric to approximate the
optimal trajectory. The intuition is that its value is large for those
candidate locations where a measurement would be highly informative.
Collecting measurements at those locations is therefore expected to greatly
reduce the error. This means that the trajectory of the UAV must pass through
those locations where the uncertainty is high.

Finding such a trajectory is however challenged by the fact that the
uncertainty metric changes as new measurements are collected and these changes
cannot be generally predicted. Thus, the trajectory of the UAV should be
computed on-the-fly as measurements are acquired. However, updating the
trajectory with the reception of every new measurement may be too costly and
may lead to erratic behavior. A more sensible alternative is to adopt a
_receding horizon_ approach that updates the trajectory every
${{t}_{\text{upd}}}$ measurements. If ${{t}_{\text{upd}}}$ is sufficiently
small, it makes sense to assume that
${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ remains approximately constant
between consecutive updates at all grid points except where a measurement is
collected, say at ${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}_{0}}$, in which case
${{u}}_{{{{g}}}_{0}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ becomes 0\. This is
because the uncertainty at a measured location is expected to be 0. From (23),
measuring at location
${\bm{x}}_{{{\tau}}}={{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ yields

$\displaystyle{{u}}({{\tilde{\bm{{r}}}}}_{{{t}}+{{t}_{\text{upd}}}})$
$\displaystyle\approx{{u}}({{\tilde{\bm{{r}}}}}_{{{t}}+{{t}_{\text{upd}}}-1})-\frac{{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}+{{t}_{\text{upd}}}-1})}{{{{G}}}}\approx\ldots$
$\displaystyle\approx{{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})-\sum_{{{\tau}}={{t}}}^{{{t}}+{{t}_{\text{upd}}}-1}\frac{{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{\tau}}})}{{{{G}}}}.$
(24)

Under this approximation, those trajectories for which
${{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ decreases quickly are those where
measurements are collected at locations ${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$
with high ${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$, which agrees with
the intuition presented at the beginning of this section. If
${{t}_{\text{upd}}}$ were so large that it were possible to collect a
measurement at all grid locations, then one could think of formulating this
problem as a discounted reward travelling salesman problem [41]. However, this
kind of problems is known to be NP-hard; see e.g. [41, 42] and references
therein. Thus, one could think of applying for example the heuristic in [41].
Nonetheless, such a task involves high complexity, which would not be
motivated given the approximation error entailed by adopting such a large
${{t}_{\text{upd}}}$; recall that the approximation in (24) applies only for
small ${{t}_{\text{upd}}}$. This would also mean that the entire trajectory
would be planned using the initial uncertainty metric values and therefore it
would not be adaptive. Besides, its complexity would render this approach
inappropriate for real-time UAV operations with limited energy and
computational power.

Thus, it is preferable to adopt an alternative approach which alternatingly
updates the uncertainty metric and plans a short trajectory through areas of
high uncertainty. To plan this trajectory, the idea here will be to set a
destination as the grid point with the highest local uncertainty and then plan
a route to that destination through points with large uncertainty. Assume
without loss of generality that the destination and trajectory are computed at
${{t}}={{t}}_{{\text{I}}}$. Thus, the destination is given by
${\bm{x}}_{{\text{F}}}\triangleq{{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}^{*}}$,
where

$\displaystyle{{{g}}}^{*}\triangleq\arg\max_{{{{g}}}\in{\bar{\mathcal{B}}}}~{}{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}_{{\text{I}}}}).$
(25)

To design the route from the current location to ${\bm{x}}_{{\text{F}}}$, it
is necessary to define a cost function that promotes trajectories through high
uncertainty locations. Let ${c}_{{t}}({\bm{x}})$ represent the cost associated
with location ${\bm{x}}$ after observing ${{\tilde{\bm{{r}}}}}_{{{t}}}$ and
defined in such a way that the higher the uncertainty value, the lower the
cost at that location. Since the uncertainty values provided by the estimators
in Secs. III and IV are associated only with grid points, it is natural to set
the cost of a location based on the uncertainty of the nearest grid point. In
particular, the cost ${c}_{{t}}({\bm{x}})$ will be set to be
${c}_{{t}}({\bm{x}})={\phi}({{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))$,
where
${{{g}}}=\arg\min_{{{{g}}}\in{\bar{\mathcal{B}}}}\|{\bm{x}}-{{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}\|_{2}$
and ${\phi}(.)$ is a non-negative decreasing function. With this definition,
one could in principle set the cost of a trajectory
${\bm{X}}=[{\bm{x}}_{{{t}}_{{\text{I}}}},\ldots,{\bm{x}}_{{{t}}_{{\text{F}}}}]$
to be
$\sum_{{{\tau}}={{t}}_{{\text{I}}}}^{{{t}}_{{\text{F}}}}{c}_{{t}}({\bm{x}}_{{{\tau}}})$,
which would assign a high cost to those trajectories through low uncertainty
areas, as desired. However, this would not account for the distance between
measurement points. Thus, it is more appropriate to let this cost be

$\displaystyle{C}({\bm{X}})=\sum_{{{\tau}}={{t}}_{{\text{I}}}+1}^{{{t}}_{{\text{F}}}}\int_{{\bm{x}}_{{{\tau}}-1}}^{{\bm{x}}_{{{\tau}}}}{c}_{{{t}}}({\bm{x}})d{\bm{x}}.$
(26)

Although this cost may be convenient for finding trajectories through regions
of high uncertainty, it may well happen that the trajectory becomes too wiggly
and, possibly, too long. This means that the values of the uncertainty metric
used when planning the trajectory may become obsolete, i.e., the approximation
error of (24) may become too large, which would lead to highly suboptimal
trajectories. For this reason, it may be sensible to account also for the time
it takes to reach the destination. To this end, let ${{T}}({\bm{X}})$ be the
time that the UAV needs to follow a trajectory ${\bm{X}}$. A reasonable
simplification is that the UAV moves at constant speed ${v}$ and, therefore,

$\displaystyle{{T}}({\bm{X}})=\frac{1}{{v}}\sum_{{{\tau}}={{t}}_{{\text{I}}}+1}^{{{t}}_{{\text{F}}}}\|{\bm{x}}_{{{\tau}}}-{\bm{x}}_{{{\tau}}-1}\|_{2}.$
(27)

Observe that if one minimized ${{T}}(\cdot)$ rather than ${C}(\cdot)$, the
resulting trajectory would be the shortest path between the current position
and the destination ${\bm{x}}_{\text{{\text{F}}}}$. Thus, there clearly exists
a trade-off between minimizing ${{T}}(\cdot)$ and ${C}(\cdot)$. A “sweet-spot”
can be found by properly choosing ${\bm{X}}$ to minimize their weighted sum

$\displaystyle\mathop{\text{minimize}}_{{\bm{X}}}~{}~{}(1-{\eta}){{T}}({\bm{X}})+{\eta}{C}({\bm{X}})$
(28) $\displaystyle\text{s.t.:
}~{}{\bm{x}}_{{{t}}_{{\text{I}}}}={\bm{x}}_{{\text{I}}},~{}{\bm{x}}_{{{t}}_{{\text{F}}}}={\bm{x}}_{{\text{F}}},~{}{\bm{x}}_{{{t}}}\notin{\mathcal{B}}^{{\mathcal{X}}},\forall{{t}},$

where ${\mathcal{B}}^{{\mathcal{X}}}$ represents the set containing the
building locations and no-fly zones, and the smaller the value of
${\eta}\in[0,1]$, the shorter the completion time but also the lower the total
uncertainty of the trajectory. Since the UAV is flying with constant velocity
${v}$, problem (28) reads as

$\displaystyle\mathop{\text{minimize}}_{{\bm{X}}}~{}\sum_{{{\tau}}={{t}}_{{\text{I}}}+1}^{{{t}}_{{\text{F}}}}\int_{{\bm{x}}_{{{\tau}}-1}}^{{\bm{x}}_{{{\tau}}}}\Big{(}\frac{1}{{v}}(1-{\eta})+{\eta}{c}_{{t}}({\bm{x}})\Big{)}d{\bm{x}}$
(29) $\displaystyle\text{s.t.:
}~{}{\bm{x}}_{{{t}}_{{\text{I}}}}={\bm{x}}_{{\text{I}}},~{}{\bm{x}}_{{{t}}_{{\text{F}}}}={\bm{x}}_{{\text{F}}},~{}{\bm{x}}_{{{t}}}\notin{\mathcal{B}}^{{\mathcal{X}}},\forall{{t}}.$

To solve (29), it is convenient to discretize the set of candidate measurement
locations ${\bm{x}}_{{{\tau}}}$, as customary in contemporary route planners
[43, 44]. Although more general sets can be used, for simplicity, it is
assumed from now on that the measurement locations must lie on the grid
${\mathcal{G}}$. In that case, if ${\bm{x}}_{{{\tau}}}$ and
${\bm{x}}_{{{\tau}}-1}$ are adjacent horizontally, vertically, or diagonally,
then the integral term in (29) is given by

$\displaystyle\int_{{\bm{x}}_{{{\tau}}-1}}^{{\bm{x}}_{{{\tau}}}}\Big{(}\frac{1}{{v}}(1-{\eta})+{\eta}{c}_{{t}}({\bm{x}})\Big{)}d{\bm{x}}=$
$\displaystyle\|{\bm{x}}_{{{\tau}}}-{\bm{x}}_{{{\tau}}-1}\|_{2}\Big{(}\frac{1}{{v}}(1-{\eta})+\frac{{\eta}}{2}({c}_{{t}}({\bm{x}}_{{{\tau}}})+{c}_{{t}}({\bm{x}}_{{{\tau}}-1}))\Big{)}$
(30)

since, by definition, ${c}_{{t}}(\cdot)$ is a piecewise constant function.

With this simplification, problem (29) becomes a shortest path problem and,
therefore, can be solved with classical algorithms such as the well-known
Bellman-Ford algorithm [32]. The transition costs are given by (V-A). To sum
up, the UAV iteratively sets a destination using (25) and plans a trajectory
using (29). Observe that it is possible that the number of measurements of the
resulting trajectory is greater than ${{t}_{\text{upd}}}$ since the opposite
is not enforced by (29). In that case, the UAV may just follow the first
${{t}_{\text{upd}}}$ points indicated by the trajectory and replan afterwards.
This would ensure that the uncertainty values do not become obsolete.

Finally, depending on the adopted estimator, it is possible that the
uncertainty metric may significantly change between updates. Since the
trajectory will depend on the uncertainty metric, it may be convenient to
smooth ${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ to avoid sudden changes
in the trajectory. To this end, one can perform a running average as follows:

$\displaystyle{{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})={{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}){\beta}+{{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}-1})(1-{\beta}),$
(31)

where ${\beta}\in(0,1]$ is a smoothing factor. Clearly, the lower ${\beta}$,
the more dependent on historical uncertainty is
${{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$.

##  VI Experiments

This section presents numerical results to evaluate the performance of the
proposed algorithms.

Spectrum surveying is carried out with a UAV that flies at a constant height
to construct 2D power maps (${{{m}}}=2$). The region of interest is
discretized into an ${{N_{y}}\times{N_{x}}}=32\times 32$ rectangular grid with
grid point spacing of $3\text{ m}$ in both dimensions. Two datasets are
considered. In the Gudmundson dataset, the receive power is generated using
the model in Sec. III-A using two transmitters with power ${P}_{\text{Tx}}=10$
dBm and height $20$ m placed uniformly at random in ${\mathcal{X}}$. The path
loss exponent is 2, the carrier frequency is $2.4$ GHz, and antennas are
isotropic. To focus on the effects of the shadowing, ${{\sigma^{2}_{{{w}}}}}$
and ${{\sigma^{2}_{{{z}}}}}$ are set to $0$. The shadowing component is
generated with ${\delta}_{{{s}}}=50$ m, ${{\sigma^{2}_{{{{s}}}}}}=10$, and
${{{\mu}}}_{{\bar{{s}}}}=0$. The other dataset, referred to as the Rosslyn
dataset, is the one from [26], where maps are generated using the commercial
3D ray-tracing software Wireless Insite. The so-called _urban canyon_ model is
used with 6 reflections and 2 diffractions. 42 large maps (40 for training and
2 for testing) are generated in a square region of $700$ m side in the
downtown of Rosslyn, Virginia, each one with a transmitter operating at $2.4$
GHz and placed at a different location. To generate each training or testing
map, an ${{N_{y}}\times{N_{x}}}$ patch is extracted from two of these large
maps at the same position and added. For both datasets, measurements off the
grid are obtained via cubic spline interpolation.

DRUE was implemented in TensorFlow with a latent dimension of $4\times 4=16$.
The Adam optimization algorithm [45] is used for training the network with a
constant learning rate of $10^{-5}$ and a batch size of 64\. Initially, DRUE
is trained with ${\alpha}=0.5$, then with ${\alpha}=0$ (i.e., only layers of
${{g}^{{\bm{R}}}_{{\bm{{\theta}}}_{{\bm{R}}}}}$ are trainable), and finally
with ${\alpha}=1$ (i.e., only layers of
${{g}^{{\bm{V}}}_{{\bm{{\theta}}}_{{\bm{V}}}}}$ are trainable). Furthermore,
for each training map, the number of measurements is drawn uniformly at random
between $1$ and $100$ with measurement locations drawn uniformly at random
over the grid.

To evaluate the performance of the considered schemes, two metrics are
utilized. One is the root mean square error (RMSE)

$\displaystyle\text{RMSE}_{{{t}}}=\sqrt{\frac{{}\mathop{\mathbb{E}}\\{\|\bm{D}_{{\bar{\mathcal{B}}}}({{\bm{{r}}}}^{{\mathcal{G}}}-{\hat{\bm{{r}}}}^{{\mathcal{G}}}_{{{t}}})\|^{2}_{2}\\}}{\mathop{\rm
Tr}(\bm{D}_{{\bar{\mathcal{B}}}})}}$ (32)

where $\bm{D}_{{\bar{\mathcal{B}}}}$ is a diagonal matrix with ones and zeros
that sets to 0 all entries of
${{\bm{{r}}}}^{{\mathcal{G}}}-{\hat{\bm{{r}}}}^{{\mathcal{G}}}_{{{t}}}$
corresponding to grid points inside buildings or no-fly zones and
${\hat{\bm{{r}}}}^{{\mathcal{G}}}_{{{t}}}$ is given in Secs. III-C and IV-A.
The expectation averages across maps and measurement locations. The other
metric is the total uncertainty defined in (23).

###  VI-A Radio Map Reconstruction

This section corroborates the ability of the proposed schemes to estimate
radio maps. To focus on estimation aspects, no route planning is considered in
this section. Instead, measurement locations are chosen uniformly at random
throughout the area. The proposed online Bayesian and data-driven DRUE
algorithms are compared here with two competitors, namely the multikernel
algorithm in [18] and K-nearest neighbors (KNN) [46]. The multikernel
algorithm is implemented with regularization parameter $10^{-7}$ and 15
Gaussian kernels with bandwidth parameter uniformly spaced between 1 and 200
m. The KNN algorithm is implemented with $K=5$ nearest neighbors.

Fig. 4 shows the RMSEt vs. the number of measurements ${{t}}$ for the
Gudmundson dataset. Batch kriging and the online Bayesian estimator know the
mean power at each location. Recall from Sec. III that the kriging algorithm
provides MMSE estimates and, therefore, the minimum RMSE, which agrees with
Fig. 4. The online Bayesian algorithm approximates the estimates of the batch
kriging algorithm by introducing the assumption that $\tilde{{r}}_{{{t}}}$ and
${{\tilde{\bm{{r}}}}}_{{{t}}-1}$ are conditionally independent given
${{\bm{{r}}}}^{{\mathcal{G}}}$; cf. Sec. III-C. Thus, in principle, the online
algorithm should perform worse than the batch version. The fact that this is
not the case shows that the error introduced by the aforementioned assumption
is negligible and therefore empirically supports the usage of the online
algorithm.

Figure 4: Comparison of the proposed estimators with the existing state-of-
the-art for the maps generated from the Gudmundson dataset. The batch kriging
and the proposed online Bayesian estimators know the mean power of the
channels while the others do not (${\delta}_{{{s}}}=50$ m,
${{\sigma^{2}_{{{{s}}}}}}=10$, and ${{{\mu}}}_{{\bar{{s}}}}=0$). Figure 5:
Comparison of the proposed estimators with the existing state-of-the-art for
the maps generated from the Rossyln dataset. All the estimators do not know
the mean power of the channels.

The performance of DRUE is slightly worse, as expected due to the fact that it
does not know the mean. Yet, it outperforms the benchmarks.

Fig. 5 depicts the RMSE vs. the number of measurements for the Rosslyn
dataset. Two effects must be observed. First, since the locations and transmit
power of the sources are unknown, the kriging and online Bayesian algorithms
assume that the map has a zero mean. As a result, they will not provide
(approximate) MMSE estimates. In fact, they are outperformed by DRUE by a wide
margin. Second, the propagation phenomena in this dataset is more complicated
than with Gudmundson maps. The fact that DRUE offers the best performance
corroborates its ability to learn from data. Finally, although the multikernel
algorithm performs similarly to the proposed online Bayesian estimator,
observe that the former is not an online algorithm and, therefore, its
computational complexity becomes eventually unaffordable if a sufficiently
large number of measurements is collected.

###  VI-B Uncertainty Learning with DRUE

If the data distribution were Gaussian, one could assess how well the
uncertainty is learned by DRUE by evaluating the log-likelihood of test data
[36]. However, this is not applicable for the Rosslyn dataset as data there is
not Gaussian; cf. Sec. IV-B. Therefore, an alternative approach is explored
here. The idea is to plot a histogram of the normalized residuals at locations
with no associated measurements. Specifically, the normalized residual
${r}_{{{{g}}}}^{\prime}$ at an unobserved location
${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ is given by

$\displaystyle{r}_{{{{g}}}}^{\prime}=\frac{[{{\bm{{r}}}}^{{\mathcal{G}}}]_{{{{g}}}}-[{\hat{\bm{{r}}}}^{{\mathcal{G}}}_{{{t}}}]_{{{{g}}}}}{[\text{vec}({{f}^{{\bm{V}}}_{{\bm{{\theta}}}}}({\check{{\bm{R}}}}_{{{t}}}))]_{{}_{{{{g}}}}}},$
(33)

where ${{t}}$ is fixed and ${{{g}}}\notin{\mathcal{B}}\cup{\Omega}$.

Figure 6: Histogram of the normalized residual with the Rosslyn dataset, where
estimates of the power map and uncertainty metric are obtained by DRUE with
$|{\Omega}|=30$ locations drawn uniformly at random.

Recall that the network is trained so that the denominator in (33) captures
the magnitude of the error, which is given by the numerator. The histogram
plot of ${r}_{{{{g}}}}^{\prime}$ is therefore expected to be centered near $1$
if the uncertainty metric is estimated well. This is seen to be the case in
Fig. 6, which shows the normalized histogram of ${r}_{{{{g}}}}^{\prime}$ for
the Rosslyn dataset, thereby confirming that the uncertainty is learned
satisfactorily by DRUE.

###  VI-C Spectrum Surveying Experiments

This section illustrates the operation of the proposed spectrum surveying
scheme and compares its performance against three other algorithms.

The UAV is allowed to move on a horizontal plane in one out of 8 directions
that differ by $45$ degrees: East, West, North, South, SouthEast, SouthWest,
NorthEast, NorthWest. To enable a fair comparison across route planners, the
UAV collects measurements every $7$ m on its trajectory instead of at the
waypoints. Throughout, DRUE is used to estimate the uncertainty metric for the
proposed minimum cost planner. The reciprocal decreasing function
${\phi}({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))=1/({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})+{\epsilon})$
is used to obtain the cost function, where $\epsilon$ is a small positive
constant, and ${\eta}$ is set to $0.75$.

Fig. 7 shows sample trajectories using the proposed uncertainty-aware minimum
cost planner with DRUE for the Gudmundson (Fig. 7a) and Rosslyn (Fig. 7b)
datasets. The first and second plots in each figure respectively correspond to
the true and estimated power maps. The third plot depicts the uncertainty
metric after measurements are collected at the locations indicated by red
markers. The trajectory indicated by white markers and the destination
indicated by a magenta marker are computed using this uncertainty metric. As
expected, it is observed that the resulting trajectory traverses points of
high uncertainty.

(a)

(b)

Figure 7: Radio map estimate with the proposed uncertainty-aware minimum cost
planner for a) the Gudmundson dataset (${\delta}_{{{s}}}=50$ m,
${{\sigma^{2}_{{{{s}}}}}}=10$, and ${{{\mu}}}_{{\bar{{s}}}}=0$); and b) the
Rossyln dataset using DRUE ( ${\eta}=0.75$, ${\beta}=0.25$, and
${\phi}({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))=1/({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})+{\epsilon})$).
White boxes in (b) represent buildings.

The rest of this section compares the proposed scheme against three other
algorithms. Since this is the first work to address _adaptive_ measurement
collection, there are no competing algorithms to compare with. Instead, the
proposed uncertainty-aware trajectory planning algorithm is compared against
the non-adaptive approach in [21] and two more benchmarks. Specifically, the
considered benchmarks are the following: i) The _grid planner_ from [21],
which traverses the grid points by moving column by column. ii) A _spiral grid
planner_ , which starts in the top-left corner and moves in a rectangular
spiral fashion towards the center of the grid. iii) An _independent uniform
planner_ , which randomly selects a grid point and proceeds towards there in a
straight line. In all cases, obstacles are avoided by the trajectory.

Fig. 8 shows RMSEt and the total uncertainty (23) obtained via Monte Carlo
simulation for the aforementioned route planners with DRUE using the Rosslyn
dataset, which serves as a proxy for real-world propagation circumstances.

Figure 8: Comparison of the proposed uncertainty-aware minimum cost route
planner with other three benchmarks for the Rossyln dataset using DRUE;
${\eta}=0.75$, ${\beta}=0.25$, ${{t}_{\text{upd}}}=7$, and
${\phi}({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))=1/({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})+{\epsilon})$.

It is seen that the devised approach reduces the uncertainty and error
significantly faster than the benchmarks. In fact, the proposed algorithm
requires less than 50% of the measurements required by the best benchmarks to
reach an RMSE of approximately 5 dB. It is also seen that the reduction in the
uncertainty metric parallels the reduction in the error, which further
supports that the metric in Sec. IV-B constitutes a reasonable uncertainty
metric. The grid planner performs worst in this experiment because, with the
small number of measurements considered, the UAV has only time to explore one
side of the map.

###  VI-D Route Planner Parameters

This section discusses two experiments intended to determine the best
configuration for the proposed route planner. The first experiment explores
the impact of the selection of the non-negative decreasing function ${\phi}$
(see Sec. V-A) on the uncertainty reduction.

Figure 9: Comparison of two different cost functions with the shortest path
cost function for the Rossyln dataset using the proposed route planner, and
DRUE; ${\eta}$ = 1, ${{t}_{\text{upd}}}=7$, and ${\beta}=0.25$.

The reciprocal function
${\phi}({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))=1/({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})+{\epsilon})$
and the exponential function
${\phi}({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))=\text{exp}(-{{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))$
are compared with the shortest path planner, which results from setting
${\phi}({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))=1$. From Fig. 9,
the steepest uncertainty and RMSE reduction are obtained with the reciprocal
function.

Fig. 10 investigates the trade-off between minimizing distance and maximizing
uncertainty collection, which is controlled by ${\eta}$; cf. Sec. V-A. When
${\eta}=0$, the UAV follows the shortest path without considering uncertainty.
From Fig. 10, this yields the worst performance. In contrast, when ${\eta}=1$,
the focus is on traversing areas with high uncertainty. This does not result
in the fastest uncertainty reduction either. The “sweet spot” is seen to be
attained for an intermediate value of ${\eta}$, namely ${\eta}=0.75$.

Figure 10: Uncertainty decrements in the map from the Rossyln dataset for
different values of ${\eta}$ with
${\phi}({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))=1/({{\bar{u}}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})+{\epsilon})$
using DRUE. The selection of ${\eta}$ trade-off between the cost and the
distance for the trajectory planning; ${{t}_{\text{upd}}}=7$, and
${\beta}=0.25$.

##  VII Related Works

This section describes the connections between the present paper and the
closely related works.

In [36], the received power at every location is modeled as an independent
Gaussian random variable and a DNN is used to learn its parameters. The
variance estimate that it provides could be used as an uncertainty metric. In
contrast, DRUE provides an uncertainty metric applicable even when the data is
not Gaussian; cf. Fig. 1. The work in [21] utilized moving sensors to collect
measurements on a path specified by a system designer to estimate a radio map.
Thus, as opposed to [21], the route planner proposed here is adaptive to
measurements and radio map estimates. The work in [4] constructs a radio map
of the outage probability using UAVs. However, the trajectory is not planned
for sensing purposes. Instead, the goal is to minimize the integral of the
outage probability. Furthermore, a fixed destination is necessarily given
rather than updated based on the measurements gathered so far.

Related work in the mobile sensing literature includes [47, 48, 49], which
addressed the problem of deciding the path of a mobile robot to maximize
information collection about a spatial field while operating under a certain
set of path constraints. In [47], a spatial field is modeled as a Gaussian
Process (GP) [50] and the measurement locations are selected based on the
estimates of the parameters rather than directly on the measurements. The
uncertainty metric proposed in [48] is a sum of a function of measurement
values and a function of measurement locations. This is not general enough to
accommodate the case where the reward is given by the reduction of uncertainty
in radio map estimation since the uncertainty is a complicated function of
both measurement values and measurement locations. The work in [49] selects
sampling locations based on whether a time-varying parameter is less than some
random number, rather than on the available measurements and estimates. Thus,
none of these schemes are directly applicable to the spectrum surveying
problem.

##  VIII Conclusions

This paper proposed a spectrum surveying scheme where autonomous UAVs collect
radio measurements to construct a power map. The task of spectrum surveying
was decomposed into two steps: uncertainty-aware radio map estimation and
trajectory planning. Two estimators with complementary benefits were
developed. The first is an online Bayesian algorithm built upon the well-known
Gudmundson shadowing model. The price to be paid for its simplicity is its
inability to capture complex propagation phenomena. To bypass this limitation,
a data-driven uncertainty mapping technique was proposed. A tailor-made DNN
architecture was designed to estimate power maps as well as the associated
uncertainty even for non-Gaussian data. Using the uncertainty provided by
these estimators, a route planning algorithm was devised to determine
trajectories that pass through the most informative locations. As future work,
one could use a variational autoencoder to learn the posterior distribution of
the radio map and develop corresponding uncertainty metrics.

## References

  * [1] D. Romero, R. Shrestha, Y. Teganya, and S. P. Chepuri,  “Aerial spectrum surveying: Radio map estimation with autonomous UAVs,”  in Proc. IEEE Int. Workshop Mach. Learn. Signal Process., Espoo, Finland, Sep. 2020.
  * [2] S. Grimoud, S. B. Jemaa, B. Sayrac, and E. Moulines,  “A REM enabled soft frequency reuse scheme,”  in Proc. IEEE Global Commun. Conf., Miami, FL, Dec. 2010, pp. 819–823.
  * [3] H. B. Yilmaz, T. Tugcu, F. Alagöz, and S. Bayhan,  “Radio environment map as enabler for practical cognitive radio networks,”  IEEE Commun. Mag., vol. 51, no. 12, pp. 162–169, Dec. 2013.
  * [4] Y. Zeng, X. Xu, S. Jin, and R. Zhang,  “Simultaneous navigation and radio mapping for cellular-connected UAV with deep reinforcement learning,”  IEEE Trans. Wireless Commun., pp. 1–1, 2021.
  * [5] J. Chen and D. Gesbert,  “Optimal positioning of flying relays for wireless networks: A LOS map approach,”  in Proc. IEEE Int. Conf. Commun., Paris, France, May 2017, pp. 1–6.
  * [6] S. Zhang, Y. Zeng, and R. Zhang,  “Cellular-enabled UAV communication: A connectivity-constrained trajectory optimization perspective,”  IEEE Trans. Commun., vol. 67, no. 3, pp. 2580–2604, 2018.
  * [7] J. Chen, U. Yatnalli, and D. Gesbert,  “Learning radio maps for UAV-aided wireless networks: A segmented regression approach,”  in Proc. IEEE Int. Conf. Commun. IEEE, 2017, pp. 1–6.
  * [8] J. Chen and D. Gesbert,  “Efficient local map search algorithms for the placement of flying relays,”  IEEE Trans. Wireless Commun., vol. 19, no. 2, pp. 1305–1319, 2019\.
  * [9] A. Alaya-Feki, S. B. Jemaa, B. Sayrac, P. Houze, and E. Moulines,  “Informed spectrum usage in cognitive radio networks: Interference cartography,”  in Proc. IEEE Int. Symp. Personal, Indoor Mobile Radio Commun., Cannes, France, Sep. 2008, pp. 1–5.
  * [10] G. Boccolini, G. Hernandez-Penaloza, and B. Beferull-Lozano,  “Wireless sensor network for spectrum cartography based on kriging interpolation,”  in Proc. IEEE Int. Symp. Personal, Indoor Mobile Radio Commun., Sydney, NSW, Nov. 2012, pp. 1565–1570.
  * [11] A. Agarwal and R. Gangopadhyay,  “Predictive spectrum occupancy probability-based spatio-temporal dynamic channel allocation map for future cognitive wireless networks,”  Trans. Emerging Telecommun. Technol., vol. 29, no. 8, pp. e3442, 2018.
  * [12] S.-J. Kim and G. B. Giannakis,  “Cognitive radio spectrum prediction using dictionary learning,”  in Proc. IEEE Global Commun. Conf., Atlanta, GA, Dec. 2013, pp. 3206 – 3211.
  * [13] D.-H. Huang, S.-H. Wu, W.-R. Wu, and P.-H. Wang,  “Cooperative radio source positioning and power map reconstruction: A sparse Bayesian learning approach,”  IEEE Trans. Veh. Technol., vol. 64, no. 6, pp. 2318–2332, Aug. 2014\.
  * [14] B. Yang, S. He, and S.-H. G. Chan,  “Updating wireless signal map with Bayesian compressive sensing,”  in Proc. ACM Int. Conf. Mod., Anal. Simu. Wireless Mobile Sys., New York, NY, USA, 2016, MSWiM ’16, pp. 310–317, Association for Computing Machinery.
  * [15] S. He and K. G. Shin,  “Steering crowdsourced signal map construction via Bayesian compressive sensing,”  in Proc. IEEE Conf. Comp. Commun. IEEE, 2018, pp. 1016–1024.
  * [16] Y. Hu, W. Zhou, Z. Wen, Y. Sun, and B. Yin,  “Efficient radio map construction based on low-rank approximation for indoor positioning,”  Math. Probl. Eng., vol. 2013, 2013.
  * [17] D. Romero, S-J. Kim, G. B. Giannakis, and R. López-Valcarce,  “Learning power spectrum maps from quantized power measurements,”  IEEE Trans. Signal Process., vol. 65, no. 10, pp. 2547–2560, May 2017.
  * [18] J.-A. Bazerque and G. B. Giannakis,  “Nonparametric basis pursuit via kernel-based learning,”  IEEE Signal Process. Mag., vol. 28, no. 30, pp. 112–125, Jul. 2013\.
  * [19] Y. Teganya, D. Romero, L. M. Lopez-Ramos, and B. Beferull-Lozano,  “Location-free spectrum cartography,”  IEEE Trans. Signal Process., vol. 67, no. 15, pp. 4013–4026, Aug. 2019.
  * [20] M. Tang, G. Ding, Q. Wu, Z. Xue, and T. A. Tsiftsis,  “A joint tensor completion and prediction scheme for multi-dimensional spectrum map construction,”  IEEE Access, vol. 4, pp. 8044–8052, Nov. 2016.
  * [21] G. Zhang, X. Fu, J. Wang, X.-L. Zhao, and M. Hong,  “Spectrum cartography via coupled block-term tensor decomposition,”  IEEE Trans. Signal Process., vol. 68, pp. 3660–3675, 2020.
  * [22] D. Lee, D. Berberidis, and G. B. Giannakis,  “Adaptive Bayesian channel gain cartography,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., Calgary, Canada, Apr. 2018, pp. 3555–3558.
  * [23] D. Lee, S.-J. Kim, and G. B. Giannakis,  “Channel gain cartography for cognitive radios leveraging low rank and sparsity,”  IEEE Trans. Wireless Commun., vol. 16, no. 9, pp. 5953–5966, 2017\.
  * [24] Y.-Q. Xu, B. Zhang, G. Ding, B. Zhao, S. Li, and D. Guo,  “Radio environment map construction based on spatial statistics and Bayesian hierarchical model,”  IEEE Trans. Cog. Commun. Net., pp. 1–1, 2021.
  * [25] B. A. Jayawickrama, E. Dutkiewicz, I. Oppermann, G. Fang, and J. Ding,  “Improved performance of spectrum cartography based on compressive sensing in cognitive radio networks,”  in Proc. IEEE Int. Commun. Conf., Budapest, Hungary, Jun. 2013, pp. 5657–5661.
  * [26] Y. Teganya and D. Romero,  “Deep completion autoencoders for radio map estimation,”  IEEE Trans. Wireless Commun., 2021.
  * [27] X. Han, L. Xue, F. Shao, and Y. Xu,  “A power spectrum maps estimation algorithm based on generative adversarial networks for underlay cognitive radio networks,”  Sensors, vol. 20, no. 1, pp. 311, Jan. 2020.
  * [28] S. Shrestha, X. Fu, and M. Hong,  “Deep generative model learning for blind spectrum cartography with NMF-based radio map disaggregation,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. IEEE, 2021, pp. 4920–4924.
  * [29] K. Saito, Y. Jin, C. Kang, J.-I. Takada, and J.-S. Leu,  “Two-step path loss prediction by artificial neural network for wireless service area planning,”  IEICE Commun. Express, vol. 1–6, Sep. 2019.
  * [30] C. E. Rasmussen and C. K. I. Williams,  Gaussian Processes for Machine Learning, vol. 2,  MIT Press, 2006.
  * [31] M. Gudmundson,  “Correlation model for shadow fading in mobile radio systems,”  Electron. Letters, vol. 27, no. 23, pp. 2145–2146, Nov. 1991.
  * [32] R. Bellman,  “On a routing problem,”  Quart. Appl. Math., vol. 16, no. 1, pp. 87–90, 1958.
  * [33] C. M. Bishop,  Pattern Recognition and Machine Learning,  Information Science and Statistics. Springer, 2006.
  * [34] T. W. Simpson, T. M. Mauery, J. J. Korte, and F. Mistree,  “Kriging models for global approximation in simulation-based multidisciplinary design optimization,”  AIAA J., vol. 39, no. 12, pp. 2233–2241, 2001.
  * [35] S. M. Kay,  Fundamentals of Statistical Signal Processing, Vol. I: Estimation Theory,  Prentice-Hall, 1993.
  * [36] E. Krijestorac, S. Hanna, and D. Cabric,  “Spatial signal strength prediction using 3D maps and deep learning,”  in Proc. IEEE Int Conf. Commun. IEEE, 2021, pp. 1–6.
  * [37] I. Goodfellow, Y. Bengio, and A. Courville,  Deep Learning,  MIT press, 2016.
  * [38] M. Ribeiro, A. E. Lazzaretti, and H. S. Lopes,  “A study of deep convolutional auto-encoders for anomaly detection in videos,”  Pattern Recognition Letters, vol. 105, pp. 13–22, Apr. 2018.
  * [39] B. Xu, N. Wang, T. Chen, and M. Li,  “Empirical evaluation of rectified activations in convolutional network,”  arXiv preprint arXiv:1505.00853, 2015.
  * [40] D. Clevert, T. Unterthiner, and S. Hochreiter,  “Fast and accurate deep network learning by exponential linear units (ELUs),”  arXiv preprint arXiv:1511.07289, 2015.
  * [41] A. Blum, S. Chawla, D. R. Karger, T. Lane, A. Meyerson, and M. Minkoff,  “Approximation algorithms for orienteering and discounted-reward TSP,”  SIAM J. Comp., vol. 37, no. 2, pp. 653–670, 2007.
  * [42] L. N. Jerome and J. P. George,  “On trajectory optimization for active sensing in Gaussian process models,”  in Proc. IEEE Conf. Dec. and Cont. IEEE, 2009, pp. 6286–6292.
  * [43] A. Bacha, C. Bauman, R. Faruque, M. Fleming, C. Terwelp, C. Reinholtz, D. Hong, A. Wicks, T. Alberi, D. Anderson, et al.,  “Odin: Team VictorTango’s entry in the DARPA urban challenge,”  J. Field Robot., vol. 25, no. 8, pp. 467–492, 2008.
  * [44] C. Urmson, J. Anhalt, D. Bagnell, C. Baker, R. Bittner, M. Clark, J. Dolan, D. Duggins, T. Galatali, C. Geyer, et al.,  “Autonomous driving in urban environments: Boss and the urban challenge,”  J. Field Robot., vol. 25, no. 8, pp. 425–466, 2008.
  * [45] D. P. Kingma and J. Ba,  “Adam: A method for stochastic optimization,”  arXiv preprint arXiv:1412.6980, 2014.
  * [46] B. Altintas and T. Serif,  “Improving RSS-based indoor positioning algorithm via k-means clustering,”  in Proc. European Wireless 2011-Sustainable Wireless Tech. VDE, 2011, pp. 1–5.
  * [47] K.-C. Ma, L. Liu, and G. S. Sukhatme,  “Informative planning and online learning with sparse Gaussian processes,”  in Proc. IEEE Int. Conf. Robot. Auto. IEEE, 2017, pp. 4292–4298.
  * [48] C. K. Ling, K. H. Low, and P. Jaillet,  “Gaussian process planning with Lipschitz continuous reward functions: Towards unifying Bayesian optimization, active learning, and beyond,”  in AAAI Conf. Artif. Intell., 2016.
  * [49] M. Popovic, G. Hitz, J. Nieto, I. Sa, R. Siegwart, and E. Galceran,  “Online informative path planning for active classification using UAVs,”  in Proc. IEEE Int. Conf. Robot. Auto. 2017, pp. 5753–5758, IEEE.
  * [50] C. E. Rasmussen,  “Gaussian processes in machine learning,”  in Proc. Summ. Sch. Mach. learn. Springer, 2003, pp. 63–71.