# MACHINE LEARNING PREDICTION OF TIME-VARYING RAYLEIGH CHANNELS

###### Abstract

Channel state information (CSI) rapidly becomes outdated in high mobility
scenarios, degrading the performance of wireless communication systems. In
these cases, time series prediction techniques can be applied to combat the
effects of outdated CSI. Recently, it has been shown that recurrent neural
networks (RNNs) exhibit outstanding performance in time series prediction
tasks. In this paper, we investigate the performance of RNN and long short
term memory (LSTM) predictors in a simple Rayleigh flat-fading channel. We
conduct numerical experiments to evaluate whether these machine- learning
(ML)-based predictors can outperform the optimal linear minimum mean square
error Wiener predictor. Our simulation results indicate that the considered
neural network predictors outperform the Wiener predictor for small
observation window lengths and are more robust under weak channel correlation
as well as in the presence of noise. Furthermore, we show that simple shallow
RNNs are sufficient to model Rayleigh channels over a wide range of Doppler
shifts.

Index Terms—  Channel prediction, machine learning, recurrent neural networks,
long-short term memory.

##  1 Introduction

For system data throughput to be maximized, several system parameters need to
be adapted to channel conditions. These parameters include, for example, the
modulation and coding scheme, transmit power and precoders. For optimal
performance, accurate CSI is required at the transmitter. Acquiring CSI
usually involves the transmission of pilot symbols to estimate the channel
conditions. However, in high mobility scenarios, the channel undergoes rapid
variation and hence the CSI estimated from the pilot symbols quickly becomes
outdated, a phenomenon known as channel aging [1].

A solution to this problem is to predict future channel states based on
knowledge of the past channel states. The Wiener predictor [2] is a well-
established prediction method, and it has been shown to compensate for the
effects of channel aging in multiple-input multiple-output (MIMO) systems [1],
[3]. Recently, RNNs and LSTMs have been proposed to solve the channel
prediction problem. This is because of their remarkable performance in time
series classification and forecasting problems [4], [5]. The authors of [6]
investigate narrowband channel prediction using RNNs. In [7] and [8], the
authors propose RNN channel predictors in the context of MIMO systems.
Comparisons of RNN prediction performance with predictors based on the Kalman
filter and autoregressive (AR) models have been made in [9], [10] and [11]
where the RNN predictor provided better prediction accuracy in some scenarios.
In [12], an LSTM predictor is used to compensate for the negative effects of
imperfect CSI and to ensure secure communications in a massive MIMO vehicular
system. The authors in [13] leverage LSTM and gated recurrent unit (GRU)
architectures to combat outdated CSI in a MIMO system and compare its
performance to that of an RNN predictor.

In this paper, we compare the channel prediction performance of the classical
Wiener filter and modern ML solutions through computer simulations. In
particular, we determine experimentally under which conditions the considered
ML-based predictors can outperform the Wiener predictor, which is optimal in
the mean square error (MSE) sense. In so doing, we seek to assess the common
claim that ML predictors generally outperform the classical model-based
solutions. We also investigate how deep a neural network needs to be for the
considered problem. Our results indicate that shallow RNNs sufficiently model
the channel over a wide range of Doppler shifts. Furthermore, we observe that
the RNN and LSTM predictors can provide more accurate prediction results than
the Wiener filter with only a few observed samples. For instance, for a
prediction length of 1 ms and observation window lengths smaller than 5
samples, the RNN predictor has a MSE gain of about 1 dB over the Wiener
predictor. The ML-based predictors also perform better under weak channel
correlation. For example, the RNN predictor exhibits a performance gain for
prediction lengths greater than 3 ms with an observation window of 5 samples.
The ML-based predictors are more robust in the presence of noise, providing a
performance gain even at low signal-to-noise ratios (SNRs). Nevertheless, the
considered methods perform roughly the same for some scenarios, e.g., larger
prediction lengths with only a few observed samples.

Notation: Vectors and matrices are written as lowercase and uppercase boldface
letters, respectively. The transpose and the Euclidean norm of a vector are
represented by $[\cdot]^{T}$ and $\|\cdot\|$, respectively. Moreover,
$\mathrm{j}$, $\mathbb{C}$, $\mathbb{R}$ and $\mathbb{E}$ represent the
imaginary unit, the set of complex numbers, the set of real numbers, and the
expectation operator, respectively.

##  2 System Model

We consider a wireless system in which a single-antenna access point (AP)
communicates with a single-antenna user equipment (UE). The system operates on
perfectly synchronized time-division duplexing (TDD). Hence, the base station
(BS) can estimate the channel samples from uplink reference signals and use
them to predict the future states of the downlink channel. The symbol period
is $T_{s}$ and the system bandwidth is $B=1/T_{s}$. Assuming a time-varying
flat-fading channel, the discrete-time representation of the baseband signal
received by the AP at symbol period $n$ is

$y[n]=h[n]x[n]+w[n],$ (1)

with $h[n]$ representing the channel coefficient, $x[n]$ the transmitted
symbol, and $w[n]$ the zero-mean circularly symmetric complex Gaussian
(ZMCSCG) noise term with variance $\sigma_{w}^{2}$, respectively.

We assume that the UE transmits orthogonal pilot symbols $x[n]$,
$\mathbb{E}\\{|x[n]|^{2}\\}=1$, and that the AP applies least-squares (LS)
estimation to obtain the uplink channel coefficients. The LS estimate can be
expressed as

$\tilde{h}[n]=y[n]x^{-1}[n]=h[n]+w[n]x^{-1}[n]=h[n]+\tilde{w}[n].$ (2)

Note that the variance of $\tilde{w}[n]=w[n]x^{-1}[n]$ is the same as that of
$w[n]$, i.e., $\sigma_{w}^{2}$, since $x[n]$ has unit power. The AP stores $L$
past noisy LS estimates in
$\tilde{\bm{h}}{\color[rgb]{0,0,0}[n]}=\big{[}\tilde{h}[n],\dots,\tilde{h}[n-L+1]\big{]}^{T}\in\mathbb{C}^{L}$
to predict the channel state $\ell$ symbol periods ahead, i.e.,
$\hat{h}[n+\ell]$. We define the SNR as

$\mathrm{SNR}=\frac{\|\bm{h}{\color[rgb]{0,0,0}[n]}\|^{2}}{\sigma_{w}^{2}},$
(3)

where
$\bm{h}{\color[rgb]{0,0,0}[n]}=\big{[}h[n],\dots,h[n-L+1]\big{]}^{T}\in\mathbb{C}^{L}$
is the noiseless channel vector over the past $L$ symbol periods.

The $L$ estimates are divided into two parts, i.e., $L_{\mathrm{train}}$
estimates are used to calculate the autocorrelation function and to train the
neural networks, and $L_{\mathrm{test}}$ estimates are used to test the
predictors. From $L_{\mathrm{test}}$, $N$ estimates, referred to as the
predictor order or the observation window length, are used to predict a future
channel sample. These details are further explained in Section 3. Figure 1
illustrates the parameters discussed so far.

###  2.1 Channel Model

Fig. 1: Illustration of channel prediction.

The sum-of-sinusoids (SOS) method of [14] is used to simulate correlated wide-
sense stationary (WSS) Rayleigh fading channels. The SOS method represents the
channel as a superposition of randomly-phased sinusoids. Hence, $h[n]$ can be
written as

$\displaystyle h[n]\triangleq
h(nT_{s})=\frac{1}{\sqrt{M}}\sum_{m=1}^{M}\left[x_{\mathrm{I},m}(nT_{s})+\mathrm{j}x_{\mathrm{Q},m}(nT_{s})\right],$
(4)

where $M$, $x_{\mathrm{I},m}(nT_{s})$ and $x_{\mathrm{Q},m}(nT_{s})$ represent
the number of sinusoids, the $m$-th in-phase component and the $m$-th
quadrature component, respectively. Further, $x_{\mathrm{I},m}(nT_{s})$ and
$x_{\mathrm{Q},m}(nT_{s})$ are given by

$\displaystyle x_{\mathrm{I},m}(nT_{s})=A_{m}\cos{[(2\pi
f_{\mathrm{D,max}}nT_{s}+\psi_{m})\cos{(\alpha_{m})}+\phi_{m}]},$ (5a)
$\displaystyle x_{\mathrm{Q},m}(nT_{s})=B_{m}\sin{[(2\pi
f_{\mathrm{D,max}}nT_{s}+\psi_{m})\cos{(\alpha_{m})}+\phi_{m}]},$ (5b)

with $A_{m}$ and $B_{m}$ representing Gaussian-distributed random attenuations
with zero mean and unit variance, $\alpha_{m}$ and $\phi_{m}$ the angle of
arrival and the phase shift for the $m$-th path, and $f_{\mathrm{D,max}}$ the
maximum Doppler shift in Hertz. To ensure that the simulator produces WSS
random variables, an additional phase term $\psi_{m}$ is required [14]. The
terms $\alpha_{m}$, $\phi_{m}$ and $\psi_{m}$ are uniformly distributed within
the interval $[-\pi,\pi)$. It is known that the SOS model (4) approaches the
well-known Jakes model as the number of sinusoids $M$ increases [14]. The
discrete-time autocorrelation function for the Jakes model is given by

$\displaystyle\varphi_{h{\color[rgb]{0,0,0}h}}[k]=\mathbb{E}\\{h[n]h^{\ast}[n+k]\\}=J_{o}(2\pi
f_{\mathrm{D,max}}|k|T_{s}),$ (6)

and the Doppler power spectrum by

$\displaystyle\Phi_{H}(f_{D})=\frac{1}{\pi
f_{\mathrm{D,max}}\sqrt{1-\left(\frac{f_{\mathrm{D}}}{f_{\mathrm{D,max}}}\right)^{2}}},$
(7)

where $J_{o}(\cdot)$ represents the zeroth-order Bessel function of the first
kind, and $k$ the lag index.

##  3 Channel Prediction

###  3.1 Wiener Predictor

The $N$-th order Wiener predictor consists of a linear filter that combines
$N$ past LS estimates to predict $\ell$ steps ahead. The estimated channel
sample is given by

$\displaystyle\hat{h}[n+\ell]=\sum_{k=1}^{N}w^{\ast}_{k,\ell}\tilde{h}[n-k+1],$
(8)

where $\\{w_{1,\ell},\ldots,w_{N,\ell}\\}$ are the predictor coefficients. The
Wiener predictor minimizes the MSE, defined as

$\mathrm{MSE}=\mathbb{E}\big{\\{}\big{|}h[n+\ell]-\hat{h}[n+\ell]\big{|}^{2}\big{\\}}.$
(9)

The predictor coefficients
$\bm{w}_{\ell}=[w_{1,\ell},w_{2,\ell},\dots,w_{N,\ell}]^{T}$ which minimize
(9) are given by

$\displaystyle\bm{w}_{\ell}=\bm{R}_{\tilde{h}\tilde{h},\ell}^{-1}\bm{r}_{\tilde{h}\tilde{h},\ell},$
(10)

in which $\bm{R}_{\tilde{h}\tilde{h},\ell}\in\mathbb{C}^{N\times N}$ denotes
the correlation matrix of the LS samples, and
$\bm{r}_{\tilde{h}\tilde{h},\ell}\in\mathbb{C}^{N}$ the correlation vector.
These statistics rely on the sample autocorrelation function
${\color[rgb]{0,0,0}\varphi_{\tilde{h}\tilde{h}}[k]\approx\varphi_{h{\color[rgb]{0,0,0}h}}[k]}$,
and they are defined as

$\displaystyle\footnotesize\bm{R}_{\tilde{h}\tilde{h},\ell}={\color[rgb]{0,0,0}\begin{bmatrix}\varphi_{\tilde{h}\tilde{h}}[0]&\varphi_{\tilde{h}\tilde{h}}[1]&\cdots&\varphi_{\tilde{h}\tilde{h}}[N-1]\\\
\varphi^{\ast}_{\tilde{h}\tilde{h}}[1]&\varphi_{\tilde{h}\tilde{h}}[0]&\cdots&\varphi_{\tilde{h}\tilde{h}}[N-2]\\\
\varphi^{\ast}_{\tilde{h}\tilde{h}}[2]&\varphi^{\ast}_{\tilde{h}\tilde{h}}[1]&\cdots&\varphi_{\tilde{h}\tilde{h}}[N-3]\\\
\vdots&\vdots&\ddots&\vdots\\\
\varphi^{\ast}_{\tilde{h}\tilde{h}}[N-1]&\varphi^{\ast}_{\tilde{h}\tilde{h}}[N-2]&\cdots&\varphi_{\tilde{h}\tilde{h}}[0]\end{bmatrix}},$
(11a) $\displaystyle\vskip 3.0pt plus 1.0pt minus
1.0pt\bm{r}_{\tilde{h}\tilde{h},\ell}=\big{[}{\color[rgb]{0,0,0}\varphi_{\tilde{h}\tilde{h}}[\ell]\quad\varphi_{\tilde{h}\tilde{h}}[\ell+1]\quad\cdots\quad\varphi_{\tilde{h}\tilde{h}}[N-1+\ell]}\big{]}^{T}.$
(11b)

We calculate the sample autocorrelation function for a lag $k$ using the
$L_{\mathrm{train}}$ past LS estimates as

${\color[rgb]{0,0,0}\varphi_{\tilde{h}\tilde{h}}[k]}={\color[rgb]{0,0,0}\frac{1}{L_{\mathrm{train}}}}\sum_{n=0}^{L_{\mathrm{train}}-1}\tilde{h}[n]\tilde{h}^{\ast}[n+k].$
(12)

Fig. 2: RNN (left) and LSTM (right) units with sigmoid ($\sigma$) and
hyperbolic tangent activation functions. Note that they both take in the state
$c_{t-1}$ from the previous unit and combine it with input $x_{t}$ to yield an
output $y_{t}$ and a new state $c_{t}$, which is passed on to the next unit.
The LSTM unit additionally takes in the previous input $y_{t-1}$.

###  3.2 RNN and LSTM Predictors

RNN and LSTM networks are a class of ML networks that contain feedback loops,
allowing them to retain past information [15]. Due to their ability to
maintain an internal state, they perform time series predictions with
remarkable accuracy. They are typically trained using gradient descent with
back-propagation, where a set of inputs is compared to the output based on a
cost function and the prediction errors are propagated back through the
network, leading to an update of weights and biases [10]. This process is
performed iteratively until the cost function converges to a minimum.

RNNs suffer from the so-called vanishing-gradient problem, whereby the
gradient information that is back-propagated becomes too small such that it
has little impact on the weights and biases. To solve this problem, LSTM units
include additional “gates” which control the amount of information passing
through the network [16]. Figure 2 illustrates the structure of RNN and LSTM
units.

####  3.2.1 Training the Neural Networks

In our channel prediction problem, the $L_{\mathrm{train}}$ complex-valued LS
estimates are divided into overlapping blocks of $N$ samples to be fed into
the neural networks during the training phase. Each complex-valued sample
block is split into its real and imaginary part and fed into the network as
${\color[rgb]{0,0,0}2N}$ input. In particular, if $h_{r}[\cdot]$ and
$h_{i}[\cdot]$ are the real part and the imaginary part of a channel sample,
the input vector comprising $N$ samples is given by

$\displaystyle\overline{\bm{h}}=\big{[}h_{r}[0],h_{i}[0],\cdots,h_{r}[N-1],h_{i}[N-1]\big{]}^{T}\in\mathbb{R}^{2N}.$
(13)

Note that we use the same notation as that of the Wiener predictor of order
$N$ to denote the number of input neurons $2N$. The neural network is then
trained based on all $L_{\mathrm{train}}$ samples.

####  3.2.2 Predicting the Channel Coefficients

In the prediction stage, we feed the neural networks with $N$ complex-valued
samples ($2N$ real-valued samples) from the $L_{\mathrm{test}}$ test samples,
$N\leq L_{\mathrm{test}}$, to predict $h[n+\ell]$. The neural network yields a
two-dimensional $\ell$-step prediction
$\big{[}\hat{h}_{r}[n+\ell],\hat{h}_{i}[n+\ell]\big{]}^{T}$ at its output. The
output is then re-synthesized into a complex-valued prediction as
$\hat{h}[n+\ell]=\hat{h}_{r}[n+\ell]+\mathrm{j}\hat{h}_{i}[n+\ell]$.

We emphasize that the considered ML-based predictors adjust their internal
weights and biases directly from the LS estimates. By contrast, the Wiener
predictor first needs to estimate the corresponding autocorrelation function,
calculate (10), and finally predict the future channel state. Therefore, the
neural networks are _data driven_ , while the Wiener predictor is based on
statistics.

##  4 Simulation Results

In this section, we present performance comparisons between the Wiener filter
and the ML-based predictors. We perform Monte Carlo (MC) simulations to
evaluate the performance of the considered predictors. A single MC trial
proceeds as follows:

  1. 1.

Generate a channel realization using the SOS method.

  2. 2.

Obtain $L$ noisy LS channel samples

  3. 3.

Design the predictors, i.e.,

     * •

Train the neural networks using $L_{\mathrm{train}}$ samples.

     * •

Estimate the autocorrelation function and calculate the Wiener predictor
coefficients using $L_{\mathrm{train}}$ samples.

  4. 4.

Predict the channel using $N$ estimates in the test set
${\color[rgb]{0,0,0}L_{\mathrm{test}}}$.

We generate the flat-fading channel using the SOS method with $M=200$. The
sampling period $T_{s}$ is set to 1 ms with a maximum Doppler spread
$f_{\mathrm{D,max}}$ of 100 Hz. The ML-based predictors are trained using the
Adam optimization algorithm [17] with a learning rate of 0.01, a hyperbolic
tangent activation function, and the MSE as the loss function. Dropout
regularization is added to reduce overfitting. We also use early stopping
[18], where the training process is halted if the loss function does not
improve beyond a certain threshold. For simplicity and to save computation
time, we limit our models to predict only one channel sample, as opposed to a
sequence of samples. The reported results were obtained by averaging $100$
independent MC trials.

First, we evaluate the performance of six RNN predictors to gain insight into
the number of required hidden layers $n_{o}$ and the number of required hidden
neurons $n_{h}$. We obtain 500 LS estimates and use 75% of them to train the
models, while the rest of the samples are used for testing. Batch training is
used with a batch size of 16 and the SNR is fixed at 10 dB. Figures 33 and 33
show the MSE performance of six models, with respect to Doppler shift, for
prediction lengths $\ell=2$ and $\ell=5$. Here, we fix the observation window
length to $N=5$. All models perform roughly the same for $\ell=2$, while for
$\ell=5$, the models with $n_{o}=1,2$ and $n_{h}=16$ perform the best overall,
especially for larger Doppler frequencies. Guided by this result, we set
$n_{o}=1$ and $n_{h}=16$ for the rest of the simulations (for both the RNN and
LSTM predictors), choosing one hidden layer instead of two to save on training
time. These results indicate that simple models are sufficient to predict
Rayleigh channels. We believe more complex models would be necessary to
characterize other channel models, e.g., the 3GPP TR 38.901 channel model
[19].

In the second set of experiments, we evaluate the MSE for varying observation
window lengths $N$, prediction lengths $\ell$ and SNRs. We obtain 1000 LS
estimates and use 75% for training the neural networks and for calculating the
statistics of the Wiener predictor and the rest for testing. Figure 44 shows
the MSE performance with respect to various observation window lengths $N$ for
a fixed SNR of 10 dB. For $\ell=1$, the ML-based predictors outperform the
Wiener predictor for $N<5$, with the RNN predictor having a gain of about 1 dB
compared to the Wiener predictor. The three predictors exhibit roughly the
same performance at $N=6$. The ML-based predictors can learn the channel
evolution even from a few observations while the Wiener predictor requires
more observations to obtain accurate channel statistics. For $\ell=3$, the
performance gain of the ML-based predictors is negligible.

(a)

(b)

Fig. 3: Model performance for various hidden layers $n_{o}$ and hidden neurons
$n_{h}$ for: (a) $\ell=2$; and (b) $\ell=5$. SNR is set to 10 dB in both
scenarios.

(a)

(b)

(c)

Fig. 4: Performance comparisons with respect to: (a) observation window length
(10 dB SNR); (b) prediction length (10 dB SNR, $N=5$); and (c) SNR ($N=5$).

Figure 44 display the prediction MSE with respect to various prediction
lengths $\ell$. The SNR is fixed at 10 dB and $N=5$. The gain provided by the
ML-based predictors increases as $\ell$ increases and, for $\ell>3$, the RNN
predictor shows a considerable gain over the Wiener predictor. This reveals
the robustness of ML-based predictors compared to the Wiener predictor even as
the correlation between channel samples grows weaker with increasing $\ell$.
The ML-based predictors perform better than the Wiener predictor over a wide
range of SNRs, as shown in Figure 44 (here, $N=5$ as before). The performance
gain for the ML-based predictors is more pronounced for $\ell=1$ and decreases
at high SNRs. Thus, the ML-based predictors are more robust in the presence of
noise and provide more accurate predictions than the Wiener filter even at low
SNRs.

##  5 Conclusions

We have compared the prediction performance of the Wiener predictor and ML-
based predictors in a Rayleigh flat-fading channel. We show that shallow RNN
models are sufficient to characterize the channel over a wide range of Doppler
shifts for the considered system model. Simulation results indicate that, for
small prediction lengths, ML-based predictors outperform the Wiener predictor
for small observation window lengths. They are also more robust under weak
channel correlation and in the presence of noise, i.e., at large prediction
lengths and low SNRs. These performance gains arise due to the non-linearities
of the neural networks and their data-driven nature, as opposed to the Wiener
predictor, which is linear and dependent on statistical knowledge of the
channel. However, for small observation window lengths with a fixed large
prediction length, as well as for high SNRs, the ML-based predictors and the
Wiener predictor perform roughly the same.

Our results indicate that ML-based channel prediction is a promising solution
that outperforms the classical Wiener filter in some scenarios. In future
work, we intend to investigate this approach in practical systems, e.g., 5G
NR, and intelligent reflective surface-assisted systems.

## References

  * [1] K. T. Truong and R. W. Heath, “Effects of channel aging in massive MIMO systems,” _Journal of Communications and Networks_ , vol. 15, no. 4, pp. 338–351, Aug. 2013.
  * [2] S. M. Kay, _Fundamentals of Statistical Signal Processing: Estimation Theory_. USA: Prentice-Hall, Inc., 1993\.
  * [3] S. Zhou and G. Giannakis, “How Accurate Channel Prediction Needs to be for Transmit-Beamforming With Adaptive Modulation Over Rayleigh MIMO Channels?” _IEEE Transactions on Wireless Communications_ , vol. 3, no. 4, pp. 1285–1294, Jul. 2004.
  * [4] F. Karim, S. Majumdar, H. Darabi, and S. Chen, “LSTM Fully Convolutional Networks for Time Series Classification,” _IEEE Access_ , vol. 6, pp. 1662–1669, 2018.
  * [5] Y. Liu, Z. Su, H. Li, and Y. Zhang, “An LSTM based classification method for time series trend forecasting,” in _Proc. 2019 14th IEEE Conference on Industrial Electronics and Applications (ICIEA)_ , Xi’an, China, Jun. 2019, pp. 402–406.
  * [6] W. Liu, L.-L. Yang, and L. Hanzo, “Recurrent Neural Network Based Narrowband Channel Prediction,” in _Proc. 2006 IEEE 63rd Vehicular Technology Conference_ , vol. 5, Melbourne, Australia, 2006, pp. 2173–2177.
  * [7] W. Jiang and H. D. Schotten, “Multi-Antenna Fading Channel Prediction Empowered by Artificial Intelligence,” in _Proc. 2018 IEEE 88th Vehicular Technology Conference (VTC-Fall)_ , Chicago, IL, USA, Aug. 2018, pp. 1–6.
  * [8] ——, “Neural Network-Based Channel Prediction and Its Performance in Multi-Antenna Systems,” in _Proc. 2018 IEEE 88th Vehicular Technology Conference (VTC-Fall)_ , Chicago, IL, USA, Aug. 2018, pp. 1–6.
  * [9] W. Jiang, M. Strufe, and H. D. Schotten, “Long-Range MIMO Channel Prediction Using Recurrent Neural Networks,” in _Proc. 2020 IEEE 17th Annual Consumer Communications & Networking Conference (CCNC)_. Las Vegas, NV, USA: IEEE, Jan. 2020, pp. 1–6.
  * [10] W. Jiang and H. D. Schotten, “Neural Network-Based Fading Channel Prediction: A Comprehensive Overview,” _IEEE Access_ , vol. 7, pp. 118 112–118 124, 2019.
  * [11] C. M. Morales and G. S. Eslava, “Linear and Non-Linear Channel Prediction Performance for a MIMO-OFDM System,” in _Proc. 2014 IEEE 5th Latin American Symposium on Circuits and Systems_ , Santiago, Chile, Feb. 2014, pp. 1–4.
  * [12] T. Peng, R. Zhang, X. Cheng, and L. Yang, “LSTM-Based Channel Prediction for Secure Massive MIMO Communications Under Imperfect CSI,” in _Proc. 2020 IEEE International Conference on Communications (ICC)_ , Dublin, Ireland, Jun. 2020, pp. 1–6.
  * [13] W. Jiang and H. D. Schotten, “Recurrent Neural Networks with Long Short-Term Memory for Fading Channel Prediction,” in _Proc. 2020 IEEE 91st Vehicular Technology Conference (VTC2020-Spring)_ , Antwerp, Belgium, May 2020, pp. 1–5.
  * [14] M. Pop and N. Beaulieu, “Limitations of sum-of-sinusoids fading channel simulators,” _IEEE Transactions on Communications_ , vol. 49, no. 4, pp. 699–708, Apr. 2001.
  * [15] Y. Yu, X. Si, C. Hu, and J. Zhang, “A review of recurrent neural networks: LSTM cells and network architectures,” _Neural computation_ , vol. 31, no. 7, pp. 1235–1270, 2019.
  * [16] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural Computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [17] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proc. 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015._
  * [18] F. Chollet _et al._ (2015) Keras. [Online]. Available: https://github.com/fchollet/keras
  * [19] 3GPP Radio Access Network Working Group, “Study on channel model for frequencies from 0.5 to 100 GHz (release 15),” 3GPP TR 38.901, Tech. Rep., 2018\.


# Aerial Spectrum Surveying:
Radio Map Estimation with Autonomous UAVs

###### Abstract

Radio maps are emerging as a popular means to endow next-generation wireless
communications with situational awareness. In particular, radio maps are
expected to play a central role in unmanned aerial vehicle (UAV)
communications since they can be used to determine interference or channel
gain at a spatial location where a UAV has not been before. Existing methods
for radio map estimation utilize measurements collected by sensors whose
locations cannot be controlled. In contrast, this paper proposes a scheme in
which a UAV collects measurements along a trajectory. This trajectory is
designed to obtain accurate estimates of the target radio map in a short time
operation. The route planning algorithm relies on a map uncertainty metric to
collect measurements at those locations where they are more informative. An
online Bayesian learning algorithm is developed to update the map estimate and
uncertainty metric every time a new measurement is collected, which enables
real-time operation.

Index Terms—  Radio maps, UAV communications, online estimation, route
planning, active learning.

##  1 Introduction

Radio maps find a myriad of applications in wireless communications, such as
network planning, interference coordination, power control, spectrum
management, resource allocation, handoff procedure design, dynamic spectrum
access, and cognitive radio; see e.g. [1, 2]. Recently, radio maps have
received great attention for autonomous UAV communications and operations; see
e.g. [3, 4]. These observations call for the development of a technology for
“surveying” a spatial region of interest to construct a radio map. The goal of
this paper is to address this task by collecting measurements with an
autonomous UAV.

Over the last few years, a significant body of literature has addressed the
estimation of radio maps from measurements acquired by spatially distributed
sensors, typically by some form of interpolation algorithm. This includes
kriging[5], compressed sensing [6, 7], dictionary learning [8], matrix [9] and
tensor completion [10], Bayesian models [11], kernel methods [12, 13, 14],
thin-plate spline regression [15], and deep learning [16, 17]. In the context
of UAV communications, radio map estimators have been proposed in [18]. All
these schemes assume that the measurement positions are given and, hence,
cannot decide where to measure next. Another related scheme is the one in
[19], which does decide the trajectory of a UAV. However, the criterion is to
minimize an outage metric and, thus, not tailored to construct a radio map.

This paper fills this gap by proposing _aerial spectrum surveying_ , whereby a
UAV autonomously collects measurements across the area of interest and
adaptively decides where to measure next so that the time required to attain a
prescribed estimation accuracy is approximately minimized.111Although the
focus is on UAVs, most of the ideas here can be extended to other mobile
robots such as terrestrial vehicles. To this end, the following challenges are
addressed: (i) Since there are infinitely many candidate measurement locations
in 3D space, the UAV needs to judiciously select an informative finite subset
of them. To this end, a Bayesian learning scheme is adopted to estimate the
radio map along with its uncertainty across space. Since adaptively planning
the trajectory requires updating this uncertainty metric as more measurements
are collected, an online learning algorithm with constant complexity per
measurement is developed. (ii) Given the aforementioned metric, the UAV needs
to plan a trajectory that prioritizes those points with a high uncertainty. To
cope with the combinatorial complexity involved in this kind of formulations,
two approximations are explored. The first relies on a receding horizon
formulation cast as a _discounted-reward_ travelling salesman problem, for
which polynomial complexity approximations exist [20]. Since this complexity
may still be unaffordable for real-time operation on board an UAV, a simpler
waypoint-search scheme based on a shortest-path subroutine and a suitably
designed spatial cost matrix is devised. This approach provides measurement
locations at a low complexity while accounting for uncertainty and experience.
The price to be paid is an increased suboptimality.

Sec. 2 addresses the contributions in (i) whereas Sec. 3 addresses those in
(ii). The proposed scheme is validated through simulations in Sec. 4.

_Notation:_ Boldface lowercase (uppercase) denote column vectors (matrices).
For a random vector $\bm{x}$, notation $\mathcal{N}(\bm{x}|\bm{\mu},\bm{C})$
or, its short-handed version $\mathcal{N}(\bm{\mu},\bm{C})$, denotes a
Gaussian distribution with mean $\bm{\mu}$ and covariance matrix $\bm{C}$.

##  2 Online Radio Map Learning

After presenting the model, this section formulates the problems of estimating
_power_ and _service maps_ as well as their associated uncertainty.

###  2.1 Radio Map Model

Let ${\mathcal{X}}\subset\mathbb{R}^{{{d}}}$ represent the geographical region
of interest, where ${{{d}}}$ is either 2 or 3, and consider a transmitter at
location ${\bm{x}}_{\text{Tx}}\in{\mathcal{X}}$. This transmitter may
correspond to a cellular base station. The location ${\bm{x}}_{\text{Tx}}$ as
well as the transmit power ${P}_{\text{Tx}}$ can be assumed known as base
stations in contemporary cellular networks share this information with the
users. A single transmitter is assumed to keep the notation simple, but
multiple transmitters can be readily accommodated. As usual, the power
received at ${\bm{x}}\in{\mathcal{X}}$ is given in logarithmic units by

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

where each term is explained next. ${{\bar{l}}}({\bm{x}})$ captures free-space
path loss and antenna gain. ${{\bar{{s}}}}({\bm{x}})$ is the shadowing loss,
which captures attenuation due to obstructions. With the usual log-normal
assumption, let
${{\bar{{s}}}}({\bm{x}})\sim\mathcal{N}({{{\mu}}}_{{\bar{{s}}}},{{\sigma^{2}_{{{{s}}}}}})$.
Following the empirical model in [21],
$\mathop{\textrm{Cov}}({{\bar{{s}}}}({\bm{x}}),{{\bar{{s}}}}({\bm{x}}^{\prime}))={{c}}(||{\bm{x}}-{\bm{x}}^{\prime}||)$,
where function ${{c}}$ is reparameterized here as
${{c}}({\delta})={{\sigma^{2}_{{{{s}}}}}}2^{-{\delta}/{\delta}_{0}}$ with
${\delta}_{0}$ the distance at which the correlation decays to 1/2. Finally,
${{w}}({\bm{x}})$ accounts for small-scale fading, caused by the
constructive/destructive interference between the signal paths arriving at
${\bm{x}}$, as well as additional unmodeled effects. As in [11],
${{w}}({\bm{x}})$ will be modeled as $\mathcal{N}(0,{{\sigma^{2}_{{{w}}}}})$.
Additionally, it is assumed independent of ${{w}}({\bm{x}}^{\prime})$ and
${{\bar{{s}}}}({\bm{x}}^{\prime\prime})$ for all
${\bm{x}}^{\prime},{\bm{x}}^{\prime\prime}\in{\mathcal{X}}$ with
${\bm{x}}^{\prime}\neq{\bm{x}}$. For clarity, 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}}}}$.
The deterministic component ${{l}}({\bm{x}})$ can be assumed known as
${{{\mu}}}_{{\bar{{s}}}}$ can be readily estimated from a set of measurements.

To estimate the radio map, a UAV equipped with a communication module capable
of measuring power and a GPS sensor collects measurements
$({\bm{x}}_{{{\tau}}},\tilde{{r}}_{{{\tau}}})$, ${{\tau}}=0,1,\ldots$, where
$\tilde{{r}}_{{{\tau}}}\triangleq{r}({\bm{x}}_{{{\tau}}})+{{z}}_{{{\tau}}}$ is
the _received signal strength_ at ${\bm{x}}_{{{\tau}}}\in{\mathcal{X}}$ and
${{z}}_{{{\tau}}}\sim\mathcal{N}(0,{{\sigma^{2}_{{{z}}}}})$ models the
measurement error, assumed independent across ${{\tau}}$ and independent of
${{w}}({\bm{x}})$ and ${{s}}({\bm{x}}^{\prime})$ for all
${\bm{x}},{\bm{x}}^{\prime}\in{\mathcal{X}}$. The measurements and their
locations up to and including time ${{t}}$ will be arranged as
${{\tilde{\bm{{r}}}}}_{{{t}}}\triangleq[\tilde{{r}}_{0},\ldots,\tilde{{r}}_{{{t}}}]^{\top}\in\mathbb{R}^{{{t}}+1}$
and
${\bm{X}}_{{{t}}}\triangleq[{\bm{x}}_{0},\ldots,{\bm{x}}_{{{t}}}]\in\mathbb{R}^{{{{d}}}\times({{t}}+1)}$.

###  2.2 Estimation Problem Formulation

This section formulates the problem of estimating two classes of fradio maps
given a collection of measurements.

#### Power Map Estimation

Given the above model, the _power map_ ${r}({\bm{x}})$ can be estimated with a
conventional Gaussian-process estimator [22, Sec. 6.4]. Unfortunately, such
_non-parametric_ approaches incur unbounded complexity as their estimates
involve the summation of one term per data point. To circumvent this effect, a
key idea here is to aggregate the information provided by all the measurements
up to and including time ${{t}}$ by the posterior of ${r}({\bm{x}})$ 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}}$.
At these points, let (cf. (2))

$\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 _batch_ version of the problem is to obtain
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}},{\bm{X}}_{{{t}}})$
given ${{\tilde{\bm{{r}}}}}_{{{t}}}$ and ${\bm{X}}_{{{t}}}$. One can then
retrieve an estimate of ${{\bm{{r}}}}^{{\mathcal{G}}}$ as the mean of this
posterior and an uncertainty metric from the covariance. However, given the
unbounded complexity that such a task may entail, it is more convenient to
address the _online_ problem of iteratively 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 $({\bm{x}}_{{{t}}},\tilde{{r}}_{{{t}}})$ with
bounded complexity per ${{t}}$.

#### Service Map Estimation

In UAV applications, rather than knowing the exact value of ${r}({\bm{x}})$,
it is often more relevant to know the set of locations ${\bm{x}}$ that the
base station can serve with a prescribed binary rate. This is necessary e.g.
to establish a command-and-control channel or to communicate application-
dependent data. Since the scheme can be readily extended to accommodate
interference, assume for simplicity that the throughput is limited by noise
and, therefore, one can regard location ${\bm{x}}$ as served if
${r}({\bm{x}})\geq{{r}}_{\min}$ for a given ${{r}}_{\min}$. Let
${{\beta}}({\bm{x}})=1$ in that case and ${{\beta}}({\bm{x}})=0$ otherwise.
The problem in this case is to find
$p({{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{{\tilde{\bm{{r}}}}}_{{{t}}},{\bm{X}}_{{{t}}})$
for each ${{{g}}}$. The online and batch versions can be phrased as before.

###  2.3 Batch and Online Bayesian Estimators

Although the focus is on online learning, the solution to the batch problem is
briefly described first to facilitate understanding. For notational
convenience, 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}$.

#### Batch Power Map Estimator

From the model embodied by (3) and (4), it can be readily shown that
${{\bm{{r}}}}^{{\mathcal{G}}}$ is conditionally independent of
${{\tilde{\bm{{r}}}}}_{{{t}}}$ given ${\bm{{{s}}}}^{{\mathcal{G}}}$. This, in
turn, implies that

$\displaystyle
p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})=\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 to lighten the notation. 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. First, the mean vectors are clearly
$\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}]=\bm{0}$ and
$\mathop{\mathbb{E}}[{{\tilde{\bm{{r}}}}}_{{{t}}}]={\bm{{{l}}}}_{{{t}}}$. For
the covariance, let
$\mathop{\textrm{Cov}}[{\bm{{{s}}}}^{{\mathcal{G}}}]\triangleq{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}$
and write
$\mathop{\textrm{Cov}}[{\bm{{{s}}}}^{{\mathcal{G}}},{{\tilde{\bm{{r}}}}}_{{{t}}}]=\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})^{\top}]=\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}(-{\bm{{{s}}}}_{{{t}}}+{\bm{{{w}}}}_{{{t}}}+\bm{{{z}}}_{{{t}}})^{\top}]=-\mathop{\mathbb{E}}[{\bm{{{s}}}}^{{\mathcal{G}}}{\bm{{{s}}}}_{{{t}}}^{\top}]\triangleq-{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}$
as well as
$\mathop{\textrm{Cov}}[{{\tilde{\bm{{r}}}}}_{{{t}}}]=\mathop{\mathbb{E}}[({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})^{\top}]=\mathop{\mathbb{E}}[(-{\bm{{{s}}}}_{{{t}}}+{\bm{{{w}}}}_{{{t}}}+\bm{{{z}}}_{{{t}}})(-{\bm{{{s}}}}_{{{t}}}+{\bm{{{w}}}}_{{{t}}}+\bm{{{z}}}_{{{t}}})^{\top}]=\mathop{\textrm{Cov}}[{\bm{{{s}}}}_{{{t}}}]+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{t}}+1}+{{\sigma^{2}_{{{z}}}}}\bm{I}_{{{t}}+1}\triangleq{{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{t}}+1}+{{\sigma^{2}_{{{z}}}}}\bm{I}_{{{t}}+1}$.
Here, the matrices ${{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}$,
${{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}$ and
${{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}}}$ can be obtained from the covariance
function ${{c}}$ introduced in Sec. 2.1. Applying [23, Th. 10.2] to this joint
distribution, it follows that
$p({\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})=\mathcal{N}({\bm{{{s}}}}^{{\mathcal{G}}}|{{\bm{\mu}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}},{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}})$,
where

$\displaystyle{{\bm{\mu}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}$
$\displaystyle=\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}}}])$
$\displaystyle=-{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}({{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{t}}+1}+{{\sigma^{2}_{{{z}}}}}\bm{I}_{{{t}}+1})^{-1}({{\tilde{\bm{{r}}}}}_{{{t}}}-{\bm{{{l}}}}_{{{t}}})$
$\displaystyle{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}$
$\displaystyle=\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}}}]$
$\displaystyle={{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}-{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}({{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{t}}+1}+{{\sigma^{2}_{{{z}}}}}\bm{I}_{{{t}}+1})^{-1}{{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}},{\bm{{{s}}}}^{{\mathcal{G}}}},$

where
${{\bm{C}}}_{{\bm{{{s}}}}_{{{t}}},{\bm{{{s}}}}^{{\mathcal{G}}}}\triangleq{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}},{\bm{{{s}}}}_{{{t}}}}^{\top}$.
Finally, applying [22, eq. (2.115)] to obtain the conditional marginal (5)
yields
$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
${{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}\triangleq{\bm{{{l}}}}_{{{t}}}-{{\bm{\mu}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}$
and
${{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}\triangleq{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{G}}}+{{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}$,
thereby solving the batch problem.

#### Online Power Map Estimator

To address the online power map estimation problem (see Sec. 2.2), it is
convenient to decompose
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})$ into
$p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1})$ and a term
that depends on the last measurement only. However, it can be easily seen that
such a factorization is not possible due to the posterior correlation among
measurements. To sidestep this difficulty, the central idea in the proposed
online learning scheme (see also Sec. 2.2) is to use ${\mathcal{G}}$ to
summarize the information of all past measurements. Mathematically, this can
be phrased 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 can be reduced by adopting a denser grid and pays off
since it enables online estimation.

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}}}),$

where $\propto$ denotes equality up to a positive factor that does not depend
on ${{\bm{{r}}}}^{{\mathcal{G}}}$. As shown earlier in this section,
$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}})$.
Since $p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1})$ is
given in the online formulation (cf. Sec. 2.2), the online learning algorithm
can use
${{\bm{\mu}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}$
and ${{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}-1}}$
to obtain $p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})$.

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 [23, 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={{l}}({\bm{x}}_{{{t}}})+({{\bm{C}}}_{{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\bm{C}}}_{{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}})$
$\displaystyle\quad\times({{\bm{C}}}_{{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\sigma^{2}_{{{w}}}}}\bm{I}_{{{G}}})^{-1}({{\bm{{r}}}}^{{\mathcal{G}}}-{\bm{{{l}}}}^{{\mathcal{G}}})\triangleq{{\bm{a}}}_{{{t}}}^{\top}{{\bm{{r}}}}^{{\mathcal{G}}}+{{{b}}}_{{{t}}}$
$\displaystyle\mathop{\textrm{Var}}[\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}}}}}+{{\sigma^{2}_{{{z}}}}}-({{\bm{C}}}_{{{s}}({\bm{x}}_{{{t}}}),{\bm{{{s}}}}^{{\mathcal{G}}}}+{{\bm{C}}}_{{{w}}({\bm{x}}_{{{t}}}),{\bm{{{w}}}}^{{\mathcal{G}}}})$
$\displaystyle\quad\times({{\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}\triangleq{{{\lambda}}}_{{{t}}},$

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, it follows from [22, eq. (2.116)] that the requested posterior is
$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}}}}$
$\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].$

The sought algorithm applies these two update equations every time a new
measurement is acquired. 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}}}$.

#### Service Map Estimation

Since the service map ${{\beta}}({\bm{x}})$ is a function of ${r}({\bm{x}})$,
it is not surprising that the algorithm from the previous section can be
readily extended to obtain service maps. To this end, apply Bayes rule and
note that ${{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})$ is
deterministically solely determined by
${r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})$ to write

$\displaystyle
p({{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{{\tilde{\bm{{r}}}}}_{{{t}}})$
$\displaystyle=\int
p({{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}),{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})d{{\bm{{r}}}}^{{\mathcal{G}}}$
$\displaystyle=\int
p({{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{{\bm{{r}}}}^{{\mathcal{G}}},{{\tilde{\bm{{r}}}}}_{{{t}}})p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})d{{\bm{{r}}}}^{{\mathcal{G}}}$
$\displaystyle=\int
p({{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}))p({{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}})d{{\bm{{r}}}}^{{\mathcal{G}}}$
$\displaystyle=\int
p({{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}))p({r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{{\tilde{\bm{{r}}}}}_{{{t}}})d{r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}).$

Noting that
$p({{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}))=1$
if ${r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})\geq{{r}}_{\min}$ and 0
otherwise, the distribution of
${{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})$ is fully characterized by

$\displaystyle{p}^{{{\beta}}}_{{{{g}}}}\triangleq{\mathbb{P}}\left[{{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})=1|{{\tilde{\bm{{r}}}}}_{{{t}}}\right]$
$\displaystyle=\int_{{{r}}_{\min}}^{\infty}p({r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{{\tilde{\bm{{r}}}}}_{{{t}}})d{r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}).$
(7)

The latter expression can be evaluated through the cumulative distribution
function of a Gaussian random variable using the mean and variance of
$p({r}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})|{{\tilde{\bm{{r}}}}}_{{{t}}})$
obtained earlier.

##  3 Adaptive Trajectory Design

Since the UAV can navigate to arbitrary locations in ${\mathcal{X}}$ to
acquire measurements, the problem becomes how to design a trajectory such that
this acquisition is performed as efficiently as possible. Since there is a
trade-off between time and estimation performance, a more formal problem
statement would be, as described later, to minimize the time required to
obtain a map estimate with a prescribed accuracy. However, quantifying
accuracy is itself a problem since the true map is not available to the UAV.
For this reason, it is necessary to develop a suitable metric that the UAV can
compute given the measurements and prior information.

###  3.1 Uncertainty Metric

The goal of this section is to design
${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})\in[0,1]$, which denotes the
uncertainty in the target (power or service) map at
${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ after observing
${{\tilde{\bm{{r}}}}}_{{{t}}}$. If the goal is to estimate a power map, it
seems reasonable to use the posterior variance. To ensure that the resulting
metric is in $[0,1]$, one may normalize by the prior variance, since the
latter constitutes an upper bound for the posterior variance. This yields

$\displaystyle{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})=\frac{[{{\bm{C}}}_{{{\bm{{r}}}}^{{\mathcal{G}}}|{{\tilde{\bm{{r}}}}}_{{{t}}}}]_{{{{g}}},{{{g}}}}}{{{\sigma^{2}_{{{{s}}}}}}+{{\sigma^{2}_{{{w}}}}}}.$
(8)

In turn, for service map estimation, note that there is little uncertainty
when ${r}({\bm{x}})$ is known to be very large or very small: the most
uncertain points are those where ${p}^{{{\beta}}}_{{{{g}}}}$ is close to
$1/2$. This is naturally quantified by the posterior entropy of
${{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})$:

$\displaystyle{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})=-{p}^{{{\beta}}}_{{{{g}}}}\log_{2}({p}^{{{\beta}}}_{{{{g}}}})-(1-{p}^{{{\beta}}}_{{{{g}}}})\log_{2}(1-{p}^{{{\beta}}}_{{{{g}}}}).$
(9)

When there are multiple transmitters, the values of the relevant metric
(either (8) or (9)) for all transmitters can be aggregated (e.g. by averaging
or taking the maximum) to obtain a single
${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ per ${{{g}}}$.

With these point-wise uncertainty metrics, one can quantify the total
uncertainty of the map after observing ${{\tilde{\bm{{r}}}}}_{{{t}}}$ by the
spatial average
${{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})\triangleq(1/{{{G}}})\sum_{{{{g}}}=0}^{{{{G}}}-1}{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$.

###  3.2 Route Planning

The UAV may use past measurements as well as prior information about the map
to decide where to measure next. Formally,
${\bm{x}}_{{{t}}+1}={{\pi}}({{\tilde{\bm{{r}}}}}_{{{t}}},{\bm{X}}_{{{t}}})$,
where function ${{\pi}}$ is the policy that needs to be designed. Informally,
one would like that ${{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ decreases as fast as
possible over time. However, the specific criterion adopted to design
${{\pi}}$ may depend on the user’s preferences. Let ${{T}}({\bm{X}}_{{{t}}})$
denote the time that the UAV needs to follow the trajectory defined by the
points in ${\bm{X}}_{{{t}}}$. A reasonable simplification is that the UAV
moves at constant speed $v$ and, therefore,
${{T}}({\bm{X}}_{{{t}}})=\sum_{{{\tau}}=1}^{{t}}\|{\bm{x}}_{{{\tau}}}-{\bm{x}}_{{{\tau}}-1}\|/v$.
One may be, for example, interested in the ${{\pi}}$ that minimizes
$\mathop{\mathbb{E}}[{{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})]$, where ${{t}}$ and
${\bm{X}}_{{{t}}}$ are such that ${{T}}({\bm{X}}_{{{t}}})$ is below a given
upper bound. Alternatively, one could minimize
$\mathop{\mathbb{E}}[{{T}}({\bm{X}}_{{{t}}})]$ subject to an upper bound on
${{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})$. Yet another possible criterion would be
to maximize the discounted reward
$\mathop{\mathbb{E}}[\sum_{{{\tau}}=1}^{{t}}{{\gamma}}^{{{T}}({\bm{X}}_{{{\tau}}})}({{u}}({{\tilde{\bm{{r}}}}}_{{{\tau}}-1})-{{u}}({{\tilde{\bm{{r}}}}}_{{{\tau}}}))]$
with ${{\gamma}}\in(0,1)$ given. Clearly, this objective promotes trajectories
with large uncertainty improvements
${{u}}({{\tilde{\bm{{r}}}}}_{{{\tau}}-1})-{{u}}({{\tilde{\bm{{r}}}}}_{{{\tau}}})$
at the beginning.

All these formulations lead to non-convex optimization problems where the
optimization variable is the function ${{\pi}}$. Thus, it is necessary to
discretize the set of candidate measurement locations, for instance by
restricting ${\bm{x}}_{{{t}}}\in{\mathcal{G}}$. Unfortunately, even in that
case, this kind of problems can be shown to be NP-hard; see e.g. [24, 20] and
references therein. Thus, one needs to resort to approximations.

Fig. 1: Trajectory (white line) followed by the autonomous UAV in a sample
spectrum surveying operation.  Fig. 2: Comparison between the proposed minimum
cost planner and three benchmarks.

In the case of power maps, note that (8) does not depend on the measurements,
but only on their location. Therefore, a (suboptimal) trajectory can be found
in an offline fashion, for example along the lines of the algorithm in [24]
and references therein. In turn, for service maps, the metric (9) does depend
on the measurements and, therefore, the trajectory should be computed on-the-
fly, as measurements are collected. However, updating the trajectory with the
reception of every new measurement may be too costly. Besides, the presence of
expectations in the aforementioned objectives renders such a task intractable.
A more sensible alternative is to update the trajectory every
${{t}_{\text{upd}}}$ measurements, assuming that
${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ remains approximately constant
between consecutive updates at all grid points except where a measurement is
collected, in which case ${{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$
becomes 0 at that point. In other words, measuring at location
${\bm{x}}_{{{t}}}={{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ yields
${{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})\approx{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}-1})-(1/{{{G}}}){{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}-1})$.

Such a receding horizon approach could be cast as an instance of the so-called
_weighted-reward traveling salesman problem_ and a solution could be
approximated by means of the algorithm in [20], which has a polynomial
complexity. For real-time UAV operations, limited by computational power, it
may be preferable to pursue alternatives with lower complexity, yet higher
suboptimality. The alternative explored here is to select, at each trajectory
update, a destination in ${\mathcal{G}}$ with highest local uncertainty. For
rectangular ${\mathcal{G}}$, if ${{u}}_{{{{g}}}}(\tilde{{r}}_{{{t}}})$ is
organized as a matrix, this destination can be found as the maximum of such a
matrix spatially filtered by a low-pass kernel. The route to reach that
destination can be sought by minimizing the line integral of
${{u}}_{{\bm{x}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ along the trajectory, where
${{u}}_{{\bm{x}}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ denotes the uncertainty at
${\bm{x}}$. This trajectory can be approximated through a shortest-path
algorithm (e.g. Bellman-Ford) with edge cost between
${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ and
${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}^{\prime}}$ given by the reciprocal of
$\int_{{{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}}^{{{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}^{\prime}}}{{u}}_{{\bm{x}}}({{\tilde{\bm{{r}}}}}_{{{t}}})d{\bm{x}}\approx\|{{{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}^{\prime}}}-{{{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}}\|({{u}}_{{{{g}}}^{\prime}}({{\tilde{\bm{{r}}}}}_{{{t}}})-{{u}}_{{{{g}}}}({{\tilde{\bm{{r}}}}}_{{{t}}}))/2$.
This clearly promotes paths through locations with high uncertainty. The
trajectory can be recomputed periodically or after reaching each destination.
Although the resulting complexity is very low, the limitation is that wiggly
trajectories, sometimes preferable [20], are penalized by this criterion.

##  4 Numerical Experiments

This section assesses the performance of the proposed algorithms by means of
simulations. To ensure reproducibility, all the code will be made available at
the authors’ websites.

For simplicity, simulations are carried out assuming that the UAV stays at a
constant height of $20$ m and, therefore, ${{{d}}}=2$. A $30\times 25$
rectangular grid is constructed over the area of interest with a separation of
$10$ m between each pair of adjacent grid points. Two transmitters of height
$10$ m are deployed at locations drawn uniformly at random over
${\mathcal{X}}$. The true map is generated by drawing
${{\bm{{r}}}}^{{\mathcal{G}}}$ from a Gaussian distribution according to (3),
where ${{l}}({\bm{x}})$ is obtained for a path loss exponent of 2, frequency
$2.4$ GHz, and isotropic transmit antennas. The transmit power is set to
${P}_{\text{Tx}}=10$ dBm for both sources. Due to space limitations, we focus
on illustrating the effect of shadowing and, thus, ${{\sigma^{2}_{{{w}}}}}$
and ${{\sigma^{2}_{{{z}}}}}$ are set to 0. The shadowing is generated with
${\delta}_{0}=50$ m, ${{\sigma^{2}_{{{{s}}}}}}=9$, and
${{{\mu}}}_{{\bar{{s}}}}=0$. To generate measurements off the grid,
${{\bm{{r}}}}^{{\mathcal{G}}}$ is interpolated using cubic splines. To
generate the service map, ${{r}}_{\min}$ is set to 5 dBm.

The route planning algorithm described at the end of Sec. 3.2 is implemented
with a $3\times 3$ kernel of all ones. A trajectory is updated only every time
the UAV reaches the destination. This update is performed through the well-
known Bellman-Ford algorithm for shortest path. The candidate waypoints lie on
${\mathcal{G}}$ and the UAV is allowed to move in one out of 8 directions that
differ 45 degrees. The uncertainty of the maps corresponding to each
transmitter is aggregated through a $\max$ operation; cf. Sec. 3.1. Since
there is no algorithm for spectrum surveying in the literature, the proposed
method is compared against three benchmarks. Each benchmark corresponds to a
different approach to plan the trajectory. The first follows parallel lines,
thus having waypoints on a rectangular grid (grid planner); the second follows
a rectangular spiral, and the third selects the next destination uniformly at
random, then moves there straight ahead. To ensure a fair comparison, all
tested approaches collect a measurement every 5 m on their trajectory. This
means that, under the assumption of constant speed, the time required to
collect ${{t}}$ measurements is the same for all approaches. Similarly, all
approaches use the proposed online estimator.

Performance is assessed in terms of total uncertainty
${{u}}({{\tilde{\bm{{r}}}}}_{{{t}}})$ and the service error rate, which is the
fraction of grid points ${{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}}$ where
${{\beta}}({{\bm{x}}}^{{\mathcal{G}}}_{{{{g}}}})$ differs from its estimate.
Fig. 1 depicts the true and estimated power map, the true and estimated
service map, as well as the service uncertainty (cf. (9)) before starting to
measure. As observed, the regions with highest uncertainty form rings around
the sources (see Sec. 3.1). The trajectory generated by the proposed route
planner is observed to target precisely these locations. Fig. 2 compares the
reduction of the error rate and uncertainty for the tested algorithms with a
Monte Carlo simulation. It is observed that the proposed scheme results in a
steeper slope, confirming that it learns the map faster than the benchmarks.

##  5 Conclusions

This paper proposes collecting radio measurements with an autonomous UAV to
construct radio power and service maps in two steps. First, an online Bayesian
learning algorithm obtains the posterior distribution of the radio map at a
set of grid locations given all past measurements. This not only provides map
estimates but also their associated uncertainty via the posterior variance.
Second, a route planner algorithm uses the relevant uncertainty metric to plan
a trajectory along areas with high uncertainty, which naturally leads to
acquire measurements at approximately the most informative locations.

## References

  * [1] 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.
  * [2] 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.
  * [3] J. Chen and D. Gesbert,  “Optimal positioning of flying relays for wireless networks: A LOS map approach,”  in Proc. IEEE Int. Conf. Commun. IEEE, 2017, pp. 1–6.
  * [4] 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, Mar. 2018.
  * [5] 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.
  * [6] 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.
  * [7] J.-A. Bazerque and G. B. Giannakis,  “Distributed spectrum sensing for cognitive radio networks by exploiting sparsity,”  IEEE Trans. Signal Process., vol. 58, no. 3, pp. 1847–1862, Mar. 2010.
  * [8] 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.
  * [9] 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, Jun. 2017.
  * [10] 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.
  * [11] 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\.
  * [12] 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.
  * [13] 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.
  * [14] D. Romero, Donghoon Lee, and G. B. Giannakis,  “Blind radio tomography,”  IEEE Trans. Signal Process., vol. 66, no. 8, pp. 2055–2069, 2018\.
  * [15] J.-A. Bazerque, G. Mateos, and G. B. Giannakis,  “Group-lasso on splines for spectrum cartography,”  IEEE Trans. Signal Process., vol. 59, no. 10, pp. 4648–4663, Oct. 2011.
  * [16] 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.
  * [17] Y. Teganya and D. Romero,  “Data-driven spectrum cartography via deep completion autoencoders,”  in IEEE Int. Conf. Commun., Jun. 2020, arXiv:1911.12810.
  * [18] 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., Paris, France, May 2017, pp. 1–6.
  * [19] Y. Zeng, X. Xu, S. Jin, and R. Zhang,  “Simultaneous navigation and radio mapping for cellular-connected UAV with deep reinforcement learning,”  available at arXiv::2003.07574, 2020.
  * [20] A. Blum, S. Chawla, D. R. Karger, T. Lane, A. Meyerson, and M. Minkoff,  “Approximation algorithms for orienteering and discounted-reward TSP,”  SIAM J. Comput., vol. 37, no. 2, pp. 653–670, Jun. 2007.
  * [21] M. Gudmundson,  “Correlation model for shadow fading in mobile radio systems,”  Electron. Letters, vol. 27, no. 23, pp. 2145–2146, 1991.
  * [22] C. M. Bishop,  Pattern Recognition and Machine Learning,  Information Science and Statistics. Springer, 2006.
  * [23] S. M. Kay,  Fundamentals of Statistical Signal Processing, Vol. I: Estimation Theory,  Prentice-Hall, 1993.
  * [24] J. Le Ny and G. J. Pappas,  “On trajectory optimization for active sensing in gaussian process models,”  in Proc. IEEE Conf. Decision Control/Chinese Control Conf., Shanghai, China, Dec. 2009, pp. 6286–6292.