# Kriging-Based 3-D Spectrum Awareness for Radio Dynamic Zones Using Aerial
Spectrum Sensors ††thanks: This work is supported in part by the NSF PAWR
award CNS-1939334 and its associated supplement for studying National Radio
Dynamic Zones (NRDZs). The authors would like to thank Wireless Research
Center for measuring antenna patterns by using an anechoic chamber. The
datasets and post-processing scripts for obtaining the results in this
manuscript are publicly accessible at [1].††thanks: S. J. Maeng, Ozgur
Ozdemir, İ. Güvenç, and Mihail L. Sichitiu are with the Department of
Electrical and Computer Engineering, North Carolina State University, Raleigh,
NC 27606 USA (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>mlsichit@ncsu.edu).

Sung Joon Maeng, Ozgur Ozdemir, Member, IEEE, İsmail Güvenç, Fellow, IEEE, and
Mihail L. Sichitiu, Member, IEEE

###### Abstract

Radio dynamic zones (RDZs) are geographical areas within which dedicated
spectrum resources are monitored and controlled to enable the development and
testing of new spectrum technologies. Real-time spectrum awareness within an
RDZ is critical for preventing interference with nearby incumbent users of the
spectrum. In this paper, we consider a 3D RDZ scenario and propose to use
unmanned aerial vehicles (UAVs) equipped with spectrum sensors to create and
maintain a 3D radio map of received signal power from different sources within
the RDZ. In particular, we introduce a 3D Kriging interpolation technique that
uses realistic 3D correlation models of the signal power extracted from
extensive measurements carried out at the NSF AERPAW platform. Using C-Band
signal measurements by a UAV at altitudes between 30 m-110 m, we first develop
realistic propagation models on air-to-ground path loss, shadowing, spatial
correlation, and semi-variogram, while taking into account the knowledge of
antenna radiation patterns and ground reflection. Subsequently, we generate a
3D radio map of a signal source within the RDZ using the Kriging interpolation
and evaluate its sensitivity to the number of measurements used and their
spatial distribution. Our results show that the proposed 3D Kriging
interpolation technique provides significantly better radio maps when compared
with an approach that assumes perfect knowledge of path loss.

###### Index Terms:

3-D spectrum awareness, AERPAW, antenna radiation pattern, I/Q samples, LTE,
Kriging interpolation, propagation modeling, RDZ, RSRP, UAV, USRP.

##  I Introduction

As the demand for advanced wireless communication services continues to grow,
efficient use of spectrum resources is becoming increasingly vital for future
wireless technologies. Therefore, the development, testing, and evaluation of
effective mechanisms to improve spectrum efficiency and sharing have become
imperative. Although there is a considerable body of literature that examines
and analyzes spectrum sharing using theoretical models and simulations, there
is a clear need to assess these approaches in real-world deployment scenarios,
taking into account realistic propagation conditions.

In this particular context, the concept of radio dynamic zones (RDZs) emerges
as a new concept [2, 3, 4], where geographical areas with dedicated spectrum
resources are effectively managed and controlled in real-time to test new
wireless innovations. This management is achieved through the sensing of
signals entering and leaving the zone [5]. RDZs serve as testing grounds for
novel spectrum sharing concepts and emerging technologies aimed at improving
spectrum efficiency within specific deployment scenarios. In RDZs, it becomes
crucial to ensure minimal or no interference to existing incumbent users of
the spectrum. Therefore, monitoring of signal leakage to passive or active
receivers outside the RDZ becomes necessary. This requires installation and
deployment of sensors within the RDZ. Monitoring scope can include both
terrestrial areas and airspace, e.g., for coexistence with unmanned aerial
vehicles (UAVs) and satellites. By monitoring and modeling the interference
levels experienced by passive receivers in these aerial scenarios, more
efficient spectrum sharing can be achieved.

The use of radio environment maps (REMs) [6] presents an effective approach
for constructing dynamic interference maps within an RDZ, which can be
generated for each location and frequency of interest. These radio maps are
generated by collecting signal power data from deployed sensors and
incorporating their corresponding location information. However, it is often
impractical to position sensors throughout the entire RDZ area. Instead,
signal power at unknown locations can be predicted using signal processing
techniques like Kriging [7], based on measurements from nearby sparsely
deployed sensors. Kriging takes advantage of the spatial correlation between
different locations to optimize the prediction of signal power. By employing
Kriging, we can efficiently interpolate and generate a radio map of signal
power using sparsely measured datasets from the sensors.

In the existing literature, several studies have focused on modeling the
spatial correlation of shadowing in received signals [8, 9], with experimental
measurements provided in [10, 11]. The application of Kriging for generating
radio maps of signal power has been validated using both simulated and real
datasets [12]. The potential of Kriging for spectrum monitoring and
interference management has been explored in [13], while [14] extends Kriging
interpolation to spectrum interpolation and analyzes it using measurement
datasets. For ground-to-UAV communications in suburban environments, path loss
and shadowing have been modeled based on measurement datasets [15, 16].
Additionally, the spatial correlation along the linear trajectory of a UAV has
been investigated [17]. In our recent works, we introduce the RDZ concept and
discuss its features and requirements [3]. Furthermore, we propose a leakage
sensing algorithm using Kriging in the two-dimensional (2D) plane of the RDZ
[18]. Notably, to the best of our knowledge, the literature does not address
the use of Kriging to obtain a three-dimensional (3D) aerial radio map based
on measurements obtained from unmanned aerial vehicles (UAVs).

In this paper, we propose to develop and use a 3D radio map to effectively
sense signal leakage from an RDZ to the receivers outside of the RDZ. We
employ a UAV as a mobile aerial sensor, collecting signal power measurements
from distinct receivers within the RDZ. The 3D interpolation of the collected
signal power is performed using the Kriging technique. The proposed method is
thoroughly analyzed and validated through a measurement campaign. The main
contributions of this paper can be summarized as follows:

  * •

Modeling 3D Radio Propagation: Considering a 3D spectrum sensing scenario, we
develop and analyze a path loss model that accounts for spatially correlated
shadowing, two-ray wireless propagation, and _measured_ antenna radiation
patterns to accurately model 3D radio propagation. We integrate 3D antenna
measurements obtained in an anechoic chamber and study improvements in model
accuracy when compared to using dipole and omnidirectional antenna patterns.

  * •

Semi-Variogram Based Kriging Interpolation: We introduce a novel method for
Kriging interpolation specifically designed for 3D spectrum monitoring. This
approach leverages a semi-variogram technique to achieve accurate and
efficient interpolation across a 3D volume using a limited set of
measurements.

  * •

Comparison with Measurement Data: We evaluate and compare the accuracy of our
proposed 3D propagation models with the measurement data collected using
software-defined radios (SDRs) at various UAV altitudes. This analysis
provides valuable insights into the performance and reliability of the
proposed approach.

The rest of this paper is organized as follows. In Section II, we present the
system model for 3-D spectrum sensing, radio propagation, and spatial
correlation in an RDZ, while in Section III, we introduce the Kriging-based
signal interpolation method for generating a 3D radio map. In Section IV, we
describe the details of our measurement campaigns for obtaining I/Q signal
samples at a UAV from an LTE-based signal source on the ground, and our
measurements in an anechoic chamber for characterizing the antenna radiation
patterns. In Section V, we analyze the effectiveness of the proposed 3D path-
loss models in predicting the received signal power at different UAV altitudes
and locations. We present numerical results on Kriging-based 3D radio map
interpolation for various scenarios in Section VI and the last section
concludes the paper.

##  II System Model

In this section, we present the models utilized for spectrum sensing within an
RDZ. Specifically, we consider a scenario where an aerial spectrum sensor
traverses the area and captures received signals from a base station (BS).
Radio propagation, correlation, and antenna radiation pattern models are also
presented.

###  II-A 3-D Spectrum Sensing with an Aerial Mobile Sensor

An RDZ should protect incumbent users outside of the zone by controlling and
managing interference signals radiating from inside the zone. The incumbent
users may include smart devices and aerial vehicles, as well as sensitive
scientific passive receivers such as satellites and ground-based radio
astronomy receivers in radio quiet zones (RQZs) [19]. Our envisioned RDZ
concept is illustrated in Fig. 1. The real-time spectrum sensing within the
boundary of the RDZs is conducted by deployed fixed / mobile ground and aerial
sensor nodes, which is an essential technique to manage dynamic spectrum
usage. The UAV moves across the RDZ space along a multi-altitude trajectory,
capturing signal data throughout.

This paper primarily focuses on the study of real-time signal sensing in the
volume of space to monitor the signal leakage from RDZs. Mobile aerial nodes,
in the form of UAVs, collect signal power data as they follow predefined
trajectories. Subsequently, the RDZ system leverages the collected dataset
from the aerial nodes to generate a radio map depicting the signal power
surrounding the RDZ space. The interpolation of this dataset facilitates the
construction of a comprehensive representation of signal power distribution.

Figure 1: Illustration of an RDZ with aerial and ground sensors and users.

###  II-B Radio Propagation Model

The location of a BS and a UAV can be represented by

$\displaystyle\mathbf{l}^{\rm bs}$ $\displaystyle=(\psi^{\rm bs},\omega^{\rm
bs},h^{\rm bs}),\;\mathbf{l}^{\rm uav}(t)=(\psi^{\rm uav},\omega^{\rm
uav},h^{\rm uav}),$ (1)

where $\psi$, $\omega$, and $h$ denote the latitude, longitude, and altitude
of the location. Note that although the location can be generally represented
by x, y, z in 3D Cartesian coordinates, we express it by latitude, longitude,
and altitude to use the information given by GPS sensors. The time-varying
location of a UAV is given by $\mathbf{l}^{\rm uav}(t)$. The horizontal
distance and the vertical distance between a BS and a UAV can be expressed as
[20]

$\displaystyle d_{\rm h}(\mathbf{l}^{\rm bs},\mathbf{l}^{\rm uav})$
$\displaystyle=\arccos\left(\sin\psi^{\rm uav}\sin\psi^{\rm bs}\right.$
$\displaystyle\left.+\cos\psi^{\rm uav}\cos\psi^{\rm bs}\cos(\omega^{\rm
bs}-\omega^{\rm uav})\right)\times A,$ (2) $\displaystyle d_{\rm
v}(\mathbf{l}^{\rm bs},\mathbf{l}^{\rm uav})$ $\displaystyle=|h^{\rm
bs}-h^{\rm uav}|,$ (3)

where $A$ is the radius of the earth ($\approx 6378137$ m). Then, the 3D
distance between a BS and a UAV is given by

$\displaystyle d_{\rm 3D}(\mathbf{l}^{\rm bs},\mathbf{l}^{\rm uav})$
$\displaystyle=\sqrt{d_{\rm h}(l^{\rm bs},l^{\rm uav})^{2}+d_{\rm v}(l^{\rm
bs},l^{\rm uav})^{2}}.$ (4)

Next, the elevation angle between a BS and a UAV can be expressed as

$\displaystyle\theta_{l}$ $\displaystyle=\tan^{-1}\left(\frac{d_{\rm
v}}{d_{\rm h}}\right).$ (5)

Figure 2: The illustration of the two-ray ground reflection model.

To develop a propagation model, we make use of a first-order approximation and
consider the rural environment in which we collect measurements. In this
scenario, we employ the two-ray ground reflection model to represent the path
loss between a BS and a UAV. This model accounts for a line-of-sight (LoS)
path as well as a strong ground reflection path, both contributing to the
received signal as the two dominant paths in an open area such as a rural
environment. The path loss characterized by the two-ray ground reflection
model can be expressed as follows [21, Chapter 2]:

$\displaystyle\mathsf{PL}_{\rm twm}(\mathbf{l}^{\rm bs},\mathbf{l}^{\rm
uav})=\left(\frac{\lambda}{4\pi}\right)^{2}\bigg{|}\underbrace{\frac{\sqrt{\mathsf{G}_{\rm
bs}(\phi_{l},\theta_{l})\mathsf{G}_{\rm uav}(\phi_{l},\theta_{l})}}{d_{\rm
3D}}}_{\text{LoS signal}}$
$\displaystyle+\underbrace{\frac{\Gamma(\theta_{r})\sqrt{\mathsf{G}_{\rm
bs}(\phi_{r},\theta_{r})\mathsf{G}_{\rm
uav}(\phi_{r},\theta_{r})}e^{-j\Delta\tau}}{r_{1}+r_{2}}}_{\text{ground
reflected signal}}\bigg{|}^{2},$ (6)

where $\mathsf{G}_{\rm bs}(\phi,\theta)$, $\mathsf{G}_{\rm uav}(\phi,\theta)$,
$\lambda$, $\phi$ denote the antenna gain of a BS, antenna gain of a UAV,
wave-length, and azimuth angle, respectively,
$\theta_{r}=\tan^{-1}\left(\frac{h^{\rm bs}+h^{\rm uav}}{d_{\rm h}}\right)$
represents ground reflection angle, and
$\Delta\tau=\frac{2\pi(r_{1}+r_{2}-d_{\rm 3D})}{\lambda}$ indicates the phase
difference between two paths. The distance and the angle parameters in the
two-ray ground reflection model are illustrated in Fig. 2. The ground
reflection coefficient with the vertically polarized signal is given by

$\displaystyle\Gamma(\theta_{r})$
$\displaystyle=\frac{\varepsilon_{0}\sin\theta_{r}-\sqrt{\varepsilon_{0}-\cos^{2}\theta_{r}}}{\varepsilon_{0}\sin\theta_{r}+\sqrt{\varepsilon_{0}-\cos^{2}\theta_{r}}},$
(7)

where $\varepsilon_{0}$ is the relative permittivity of the ground and the
value depends on the type of the ground. Two signal components in (II-B) are
received and combined with a phase difference. If we only consider the first
LoS term in the path loss, we can obtain the free-space path loss model, given
as

$\displaystyle\mathsf{PL}_{\rm
fs}=\left(\frac{\lambda}{4\pi}\right)^{2}\left|\frac{\sqrt{\mathsf{G}_{\rm
bs}(\theta_{l})\mathsf{G}_{\rm uav}(\theta_{l})}}{d_{\rm 3D}}\right|^{2}.$ (8)

Using (II-B), the received signal power of a UAV in dB scale can be expressed
as

$\displaystyle r$ $\displaystyle=\mathsf{P}_{\rm Tx}-\mathsf{PL}_{\rm
twm}^{(\rm dB)}+w,$ (9)

where $\mathsf{P}_{\rm Tx}$, $w$ denote transmit power and shadowing
component, respectively. Note that the path loss term in (9) is converted to
dB scale. The shadowing term generally follows a lognormal distribution and is
modeled by a zero-mean Gaussian process with a spatial covariance [8]. The
correlation between received signals at two different locations is generally
characterized by the function of the distance between those locations. Note
that we do not take into account small-scale fading in the received signal
since we assume that the effect is eliminated by averaging the samples within
the proper time interval [10].

###  II-C Spatial Correlation Model of Received Signal

In this section, we focus on describing the correlation function between the
received signals at different locations of a UAV. Since the spatial
correlation primarily depends on the shadowing component ($w$) in the received
signal in (9), we can capture the correlation between received signals ($r$)
using the correlation between the shadowing components without loss of
generality. It is well-known that the correlation between two different
locations is characterized by a function of their physical distance.
Typically, this correlation exponentially attenuates as the physical distance
between the locations increases [10]. However, most existing works in the
literature primarily focus on terrestrial networks and do not fully consider
3D topologies. Due to this limitation, the spatial correlation between two
locations with different vertical positions (heights) has not been extensively
studied to our best knowledge. Considering the unique characteristics of UAV-
based scenarios, where altitude plays a crucial role, it becomes essential to
investigate and understand the spatial correlation between locations at
different vertical positions. This exploration will allow for a more
comprehensive modeling of the correlation in 3D scenarios, considering the
impact of vertical distance in addition to horizontal distance.

In our work, we first model the spatial correlation as a function of the
vertical distance ($d_{\rm v}$) as well as the horizontal distance ($d_{\rm
h}$). Then, we define the correlation function between 3D locations as a
function of both the vertical distance and the horizontal distance. The
spatial correlation between two different locations of a UAV, i.e., between
$l^{\rm uav}_{i}$ and $l^{\rm uav}_{j}$, can be expressed as

$\displaystyle R(l^{\rm uav}_{i},l^{\rm uav}_{j})=R(d_{\rm v},d_{\rm
h})=\frac{\mathbb{E}\left[w(l^{\rm uav}_{i})w(l^{\rm
uav}_{j})\right]}{\sigma_{w}^{2}},$ (10)

where $\sigma_{w}^{2}$ is the variance of shadowing. Once again, the proposed
correlation is the function of both the vertical distance and the horizontal
distance.

###  II-D Antenna Radiation Model

The antenna gain effect of a transmitter and a receiver in the received signal
is captured in the path loss model in (II-B), using $\mathsf{G}_{\rm
bs}(\phi,\theta)$, $\mathsf{G}_{\rm uav}(\phi,\theta)$. In typical terrestrial
communications, the antenna gain is simply modeled by a constant gain. This is
due to the fact that a dipole antenna is usually characterized as an omni-
directional antenna radiation pattern in the azimuth angle domain, or sectored
directional antennas make the antenna pattern mostly uniform in the azimuth
angle domain. However, air-to-ground communications require considering the
variation of the antenna gain in the elevation angle domain. The antenna
pattern in the elevation domain is typically far from being uniform and
therefore we should consider the elevation angle-dependent radiation pattern
in modeling the antenna gain.

##  III 3D Radio Map Interpolation using Kriging

In this section, we introduce an efficient radio map interpolation technique
using Kriging [13]. This method utilizes measurement data obtained from
sparsely deployed spectrum sensors within an RDZ. The interpolation process
allows us to estimate signal values at unsampled locations based on the
available measurements. We first introduce how to calculate a semi-variogram,
and subsequently, introduce our Kriging based interpolation approach for 3D
RDZ scenarios. Different than the existing Kriging techniques in the
literature, we consider the 3D geometry in spatial correlation with a portable
aerial sensor, which enables us to interpolate the radio map in a 3D volume.

###  III-A Semi-variogram

In geostatistics, the semi-variogram represents the degree of spatial
dependency on different locations which is utilized in Kriging interpolation.
The semi-variogram between a UAV’s locations $l^{\rm uav}_{i}$, $l^{\rm
uav}_{j}$ is defined as

$\displaystyle\gamma(l^{\rm uav}_{i},l^{\rm
uav}_{j})=\frac{1}{2}\text{var}\left(r(l^{\rm uav}_{i})-r(l^{\rm
uav}_{j})\right).$ (11)

If the covariance function of a stationary process exists, we can obtain the
semi-variogram from the spatial correlation in (10) as follows for our
considered 3D RDZ scenario [22]:

$\displaystyle\gamma(l^{\rm uav}_{i},l^{\rm uav}_{j})$
$\displaystyle=\frac{\sigma_{w}^{2}}{2}\left(R(l^{\rm uav}_{i},l^{\rm
uav}_{i})+R(l^{\rm uav}_{j},l^{\rm uav}_{j})-2R(l^{\rm uav}_{i},l^{\rm
uav}_{j})\right)$ $\displaystyle=\sigma_{w}^{2}\left(1-R(l^{\rm
uav}_{i},l^{\rm uav}_{j})\right)=\sigma_{w}^{2}\left(1-R(d_{\rm v},d_{\rm
h})\right),$ (12)

where $\sigma_{w}^{2}$ captures the variance of the shadowing term $w$ in (9)
as defined earlier, and $R(l^{\rm uav}_{i},l^{\rm uav}_{i})$ is as defined in
(10). We assume that $\sigma_{w}^{2}$ is constant at given set of locations
while deriving (III-A).

###  III-B Kriging Interpolation

The ordinary Kriging is the optimal prediction method in squared-error loss
from the observed data at known spatial locations where the error of the
spatial prediction of an unknown location is minimized [22]. It interpolates
the signal strength of the arbitrary locations by using the linear combination
of the signal strength of the nearby locations. The ordinary Kriging problem
can be formulated as follows [13]:

$\displaystyle\min_{\mu_{1},\dots,\mu_{M}}$
$\displaystyle\quad\mathbb{E}\left[\left(\hat{r}(\mathbf{l}^{\rm
uav}_{0})-r(\mathbf{l}^{\rm uav}_{0})\right)^{2}\right],$ (13) s.t.
$\displaystyle\quad\hat{r}(\mathbf{l}^{\rm
uav}_{0})=\sum_{i=1}^{M}\mu_{i}r(\mathbf{l}^{\rm uav}_{i}),$ (14a)
$\displaystyle\quad\sum_{i=1}^{M}\mu_{i}=1,$ (15a)

where $l^{\rm uav}_{0}$ is a location to predict an unknown parameter,
$\mu_{i}~{}(i=1,\cdots,M)$ are weighting parameters and $M$ indicates the
number of nearby measured samples to use.

The above problem can be solved by following steps [13]. First, we convert the
original problem to an equivalent Lagrange expression:

$\displaystyle\min_{\mu_{1},\dots,\mu_{M}}$
$\displaystyle\mathbb{E}\left[\left(r(l^{\rm
uav}_{0})-\sum_{i=1}^{M}\mu_{i}r(l^{\rm
uav}_{i})\right)^{2}\right]-\kappa\left(\sum_{i=1}^{M}\mu_{i}-1\right),$ (16)

where $\kappa$ denotes the Lagrange multiplier. After a few mathematical
steps, the objective function in (16) can be reformulated as

$\displaystyle\sigma^{2}_{w}+2\sum_{i=1}^{M}\mu_{i}\gamma(l^{\rm
uav}_{0},l^{\rm
uav}_{i})-\sum_{i=1}^{M}\sum_{j=1}^{M}\mu_{i}\mu_{j}\gamma(l^{\rm
uav}_{i},l^{\rm uav}_{j})$
$\displaystyle-\kappa\left(\sum_{i=1}^{M}\mu_{i}-1\right),$ (17)

where $\gamma(l^{\rm uav}_{i},l^{\rm uav}_{j})$ is as defined in (11).
Finally, we can find the optimal solution that minimizes the objective
function by the first derivative of (III-B) with respect to
$\mu_{1},\dots,\mu_{M}$, which is given by

$\displaystyle\sum_{j=1}^{M}\mu_{j}\gamma(l^{\rm uav}_{i},l^{\rm
uav}_{j})-\gamma(l^{\rm uav}_{0},l^{\rm uav}_{i})+\kappa^{\prime}=0.$ (18)

We can also express (18) as a linear matrix equation as:

$\displaystyle\left[\begin{array}[]{c c c c}\gamma(l^{\rm uav}_{1},l^{\rm
uav}_{1})&\cdots&\gamma(l^{\rm uav}_{1},l^{\rm uav}_{M})&1\\\ \gamma(l^{\rm
uav}_{2},l^{\rm uav}_{1})&\cdots&\gamma(l^{\rm uav}_{2},l^{\rm uav}_{M})&1\\\
\vdots&\vdots&\vdots&\vdots\\\ \gamma(l^{\rm uav}_{M},l^{\rm
uav}_{1})&\cdots&\gamma(l^{\rm uav}_{M},l^{\rm uav}_{M})&1\\\ 1&\cdots&1&0\\\
\end{array}\right]\left[\begin{array}[]{c}\mu_{1}\\\ \mu_{2}\\\ \vdots\\\
\mu_{M}\\\ \kappa^{\prime}\end{array}\right]$ (29)
$\displaystyle=\left[\begin{array}[]{c}\gamma(l^{\rm uav}_{0},l^{\rm
uav}_{1})\\\ \gamma(l^{\rm uav}_{0},l^{\rm uav}_{2})\\\ \vdots\\\
\gamma(l^{\rm uav}_{0},l^{\rm uav}_{M})\\\ 1\end{array}\right].$ (35)

Then, we can easily obtain the optimal $\mu_{1}^{\star},\dots,\mu_{M}^{\star}$
from (29) and interpolate the received signal powers of unknown location
$l^{\rm uav}_{0}$ by

$\displaystyle\hat{r}(l^{\rm uav}_{0})=\sum_{i=1}^{M}\mu_{i}^{\star}r(l^{\rm
uav}_{i}).$ (36)

Note that accurate characterization of the 3D semi-variogram in (11) is
critical for the interpolation in (36). The next section describes our
measurements that will be used to characterize the 3D semi-variogram.

(a) Experiment site

(b) BS tower

(c) Drone

Figure 3: (a) The area where air-to-ground propagation data has been collected
in AERPAW, (b) The fixed node tower (30 feet high) that includes the SDR
serving as the LTE eNB, and (c) the drone that carries the receiver SDR.

##  IV Measurement Campaign Overview

In this section, we describe the details of our radio propagation
measurements. We present our measurement setup, define UAV trajectory used,
and describe our approach for characterizing antenna effects.

###  IV-A Measurement Setup

The measurement campaign was conducted at the Lake Wheeler Road Field Labs
(LWRFL) site in Raleigh, NC, USA, which is one of the two sites in the NSF
Aerial Experimentation and Research Platform for Advanced Wireless (AERPAW).
The experimental area, depicted in Fig. 3a, can be classified as an open rural
environment, ensuring LoS conditions between a UAV and the BS throughout the
entire duration of the experiments. Fig. 3b and Fig. 3c present photos of the
base station (BS) tower and the drone used during the measurement campaign.
The BS tower stands at a height of 10 meters and is equipped with a single
dipole transmit antenna. On the other hand, the drone is equipped with a
vertically oriented single dipole receiver antenna and a GPS receiver to
accurately track its position. To facilitate the measurements, the srsRAN
open-source Software Defined Radio (SDR) software was utilized to implement an
LTE evolved NodeB (eNB) at the BS tower, as shown in Fig. 3b. The eNB
continuously transmitted common reference symbols (CRSs) during the
measurement campaign.

During the measurement campaign, the drone collects raw I/Q data samples using
a Software Defined Radio (SDR) that is attached to it. Specifically, the USRP
B205mini from National Instruments (NI) is utilized as the SDR device, both at
the BS tower and on the UAV. For post-processing the raw I/Q data, we employ
Matlab’s LTE toolbox. Within this toolbox, we calculate the Reference Signal
Received Power (RSRP) for each location of the UAV. To ensure efficient
processing and analysis, we collect 20 ms segments of data out of every 100
ms. Within each 20 ms segment, we extract a 10 ms duration for subsequent
post-processing. Throughout the paper, the terms “received signal” and “RSRP”
are used interchangeably to refer to the measured signal strength. The major
specifications of the transmitter and the receiver are listed in Table I.

TABLE I: Measurement Setup for Experiments BS Tower (Transmitter) |
---|---
Technology | LTE
Tower height | 10 m
Transmit power | 10 dBm
Carrier frequency | 3.51 GHz
Bandwidth | 1.4 MHz
Antenna | Dipole antenna (RM-WB1)
UAV (Receiver) |
Antenna | Dipole antenna (SA-1400-5900)
UAV heights | {30, 50, 70, 90, 110} m

###  IV-B UAV Trajectory

(a) Top view of trajectory and RSRPs ($h=110$ m). (b) 3D view of trajectory
and RSRP.

Figure 4: The trajectory of the UAV and measured RSRP in RDZ experiments.
Multiple flights with different heights are plotted in the bottom figure.

We conduct the experiments multiple times by changing the altitude (height) of
the UAV from 30 m to 110 m at increments of 20 m. In each flight, the UAV
flies an identical predefined trajectory with a different fixed height. In
particular, the UAV flies on a zig-zag pattern through the experiment site,
between south and north waypoints, and it eventually flies back to the
starting point. The top view (at $h=110$ m) and the 3D view of the UAV
trajectories along with measured RSRPs are illustrated in Fig. 4 for flight
trajectories at $30$ m, $50$ m, $70$ m, $90$ m, and $110$ m.

###  IV-C Antenna Radiation Pattern Characterization

(a) Anechoic chamber setup for Tx antenna pattern measurement. The center
point of the chamber is adjusted to the tip of the antenna by the crossed
laser lines.

(b) Measured Tx antenna pattern in Cartesian coordinates at 3.5 GHz (by the
gain in linear scale).

(c) Rx antenna pattern for the elevation angle domain at 2.4 GHz (gain in dB
scale from the specification sheet).

(d) Combined (Tx $+$ Rx) antenna pattern in 2D angle domain in dB scale. Only
the angle space in the black rectangle is used for LoS between the BS tower
and the UAV in the experiment ($\theta_{\rm l}$), and the red rectangle is
used for the ground reflection angles ($\theta_{\rm r}$).

Figure 5: The anechoic chamber setup for the 3D antenna pattern measurement
and the Tx, Rx, and combined antenna patterns we used for the analysis.

(a) RSRP vs. time (measured antenna pattern).

(b) RSRP vs. 3D distance (measured antenna pattern).

(c) RSRP vs. elevation angle (measured antenna pattern). (d) RSRP vs. time
(dipole antenna pattern).

(e) RSRP vs. 3D distance (dipole antenna pattern). (f) RSRP vs. elevation
angle (dipole antenna pattern).

(g) RSRP vs. time (onmi-directional antenna pattern). (h) RSRP vs. 3D distance
(onmi-directional antenna pattern).

(i) RSRP vs. elevation angle (onmi-directional antenna pattern).

Figure 6: RSRP fitting with different path loss models and antenna patterns in
time, distance, and elevation angle domains.

The dipole antenna used in our experiments generally exhibits omni-directional
radiation patterns in the azimuth angle domain, but oval-shaped radiation
patterns in the elevation angle domain. The radiation pattern also varies with
the carrier frequency. We obtained the antenna pattern specifications for the
Rx dipole antenna (SA-1400-5900) from the vendor’s specification sheet, and it
shows a typical donut-shaped dipole pattern that remains consistent across
different carrier frequencies [23]. Specifically, in the specification sheet,
the antenna patterns for 1.4, 1.7, 2.4, 4.4, and 5.8 GHz frequencies are
provided and all of them have similar dipole patterns. Therefore, we adopted
the 2.4 GHz frequency antenna pattern from the specification sheet for our
analysis. However, the Tx dipole antenna (RM-WB1-DN) exhibited different
elevation angle domain patterns depending on the carrier frequency and had an
asymmetric pattern that did not guarantee omni-directionality in the azimuth
angle domain [24]. Furthermore, the specification sheet did not provide the
radiation pattern for the specific carrier frequency (3.51 GHz) used in our
experiments. To obtain the exact antenna radiation pattern for the 3.51 GHz
frequency, we conducted separate measurements of the 3D antenna pattern using
an anechoic chamber facility located at wireless research center (WRC), Wake
Forest, NC.

Fig. 5a shows a photo of the setup in the anechoic chamber during the
measurement of the Tx antenna’s 3D pattern. Fig. 5b displays the output of the
antenna measurement, visualizing the antenna pattern in 3D Cartesian
coordinates. It can be observed that the antenna pattern is not purely omni-
directional in the azimuth angle domain, and the directivity in the elevation
angle domain is not straightforward. In contrast, Fig. 5c shows the elevation
angle domain antenna pattern of the Rx antenna as provided in the
specification sheet, where the antenna pattern is specified as omni-
directional with uniform gain in the azimuth domain. Fig. 5d illustrates the
combined antenna gain from the Tx and Rx antenna patterns from Fig. 5b and
Fig. 5c, respectively, represented in the azimuth and elevation angle domain.
For all UAV heights in our experiments, the LoS angles between the Tx tower
and the UAV were within the angle space covered by the black rectangular area,
while the ground reflection angles between Tx tower and the UAV were covered
by the red rectangular area, which are illustrated in Fig. 2. This implies
that the antenna pattern used for the analysis is limited to the angles within
this space.

##  V Air-to-ground Propagation Modeling and Analysis

In this section, we review how we post-process the data for correcting errors
in altitude reported by the UAV’s GPS. Subsequently, we model the measured
RSRP using different 3D propagation models that take into account two-ray
multipath model and 3D antenna pattern.

(a) CDF of RSRP.

(b) CDF of fitting error between measured RSRP and path loss models.

Figure 7: CDF of RSRP and the fitting error with path loss models by different
antenna radiation patterns. The measured antenna pattern in Fig. 5 achieves
the closest fitting to the RSRP measurements.

(a) Fitting error vs. time

(b) Fitting error vs. 3D distance

(c) Fitting error vs. elevation angle

Figure 8: Relative fitting error between measured RSRP and the two-ray path
loss model with different antenna patterns in time, distance, and elevation
angle domains.

###  V-A Post-measurement Correction of Altitude and RSRP

During the measurements, we encountered calibration errors caused by
limitations in the SDR hardware. Specifically, the Universal Software Radio
Peripheral (USRP) mounted on the UAV exhibited a power level calibration
error, resulting in a constant offset power throughout the experiment. To
address this issue, we conducted a separate experiment to measure and
determine the offset at the USRP, which was found to be 98 dB. Subsequently,
we added this offset to the calculated RSRP values obtained from subsequent
experiments, effectively compensating for the calibration offset.

Additionally, the GPS receiver carried by the UAV exhibited an altitude
mismatch. We observed an altitude drift of approximately 6 m after the UAV
landed, when compared with the initial altitude of the UAV. To rectify this
mismatch, we applied a linear compensation approach (see [25, Fig. 6]). This
involved adjusting the altitude measurements such that the altitude at the end
of the flight matched the altitude of the initial measurement. By applying
this compensation, we aimed to ensure accurate altitude data throughout the
experiment.

(a)  (b)  (c)  (d)

Figure 9: Curve fitting for the RSRP measurements at different altitudes,
considering free-space path loss and two-ray (ground-reflection) model.

(a) Fitting error vs. 3D distance (measured antenna pattern)

(b) Fitting error vs. 3D distance (dipole antenna pattern)

(c) Fitting error vs. 3D distance (constant antenna pattern)

Figure 10: Relative fitting error between measured RSRP and the two-ray path
loss model with different antenna patterns and heights in distance domains.

###  V-B Antenna Radiation Pattern Effect in Path Loss Analysis

In this subsection, we analyze the effect of antenna radiation patterns on the
path loss fitting to the RSRP from the experiments. We consider three
different antenna pattern setups for comparison: 1) Tx and Rx 3D antenna
patterns described in Section IV-C and Fig. 5; 2) the donut shape dipole
antenna pattern using the formulation for both Tx and Rx antennas; and 3)
constant azimuth and elevation antenna gain for both Tx and Rx antennas. The
dipole antenna pattern formula in the second case is given by [26]

$\displaystyle\mathsf{G}_{\rm bs}(\theta)$ $\displaystyle=\mathsf{G}_{\rm
uav}(\theta)=\frac{\cos\left(\frac{\pi}{2}\cos\theta\right)}{\sin\theta}.$
(37)

Fig. 6 and Fig. 7 provide a comprehensive analysis of the RSRP fitting results
using different antenna patterns and path loss models in (II-B), (8). In Fig.
6, the RSRP curves for a UAV height of 70 m are presented, along with the
fitting results obtained from the free space and two-ray path loss models with
different antenna patterns. It is observed that the antenna pattern described
in Section IV-C provides the best fit to the RSRP curves, while the dipole
pattern in (37) results in the worst fit. Additionally, Fig 6a highlights that
the two-ray path loss model performs better than the free space path loss
model in capturing the deep fading of RSRP.

To further evaluate the performance, Fig. 7 presents the cumulative
distribution function (CDF) of the RSRP for the 70 m height measurement, along
with the fitting results obtained from the path loss models and different
antenna patterns. The CDF of the two-ray path loss model with the antenna
pattern in Section IV-C matches closest with the CDF of the measured RSRP,
indicating a better fit. Fig. 7b shows the fitting error, which is calculated
by subtracting the measured RSRP from the fitted RSRP using the path loss
models. It is observed that the fitting error is the smallest when using the
two-ray path loss model with the antenna pattern in Section IV-C.

Fig. 8 also evaluates the fitting error with different antenna patterns in
time, distance, and elevation domains. It is observed that the dipole antenna
pattern has the largest fitting error in short and long distances. We also
observe that the fitting error is relatively high in small elevation angles.
It implies that the effect of scattering from the objects around the test site
increases the variance of the error when the elevation angle is low. Overall,
these results demonstrate that the choice of antenna pattern and path loss
model significantly impacts the accuracy of RSRP fitting for air-to-ground
communication links. The 3D antenna radiation pattern described in Section
IV-C, combined with the two-ray path loss model, provides the best fit to the
measured RSRP and minimizes the fitting error.

###  V-C Path Loss Model Fitting with Measurement

Fig. 9 illustrates the measured and fitted RSRP values as a function of 3D
distance for different UAV heights ranging from 30 m to 110 m. We adopt
measured antenna patterns in Section IV-C. The fitted curves follow the
measured RSRP values reasonably closely. It is worth noting that the two-ray
path loss model performs better in capturing the fluctuation of signal
strength due to the ground reflected path compared to the free-space path loss
model, especially when the UAV height is low. In the logarithmic scale of the
distance domain, the RSRP is expected to decrease linearly. However, in the
short distance range, a concave curve can be observed. This phenomenon is a
result of the elevation-dependent antenna gain and the dramatic change in the
elevation angle at short distances and high UAV altitudes. The 3D antenna
pattern considered in the path loss models effectively captures this effect,
leading to more accurate RSRP fitting. Overall, the results in Fig. 9
highlight the importance of considering the elevation-dependent antenna gain
and the 3D antenna pattern in accurately modeling and fitting RSRP
measurements in air-to-ground communications.

Fig. 10 shows the relative fitting error in the distance domain for all
heights with different antenna patterns. The error by the dipole antenna
pattern is relatively higher than other antenna patterns, especially when the
distance is around 100 m to 200 m due to the antenna pattern mismatch. We also
observe that the fitting error for the omnidirectional antenna pattern is
higher than the measured antenna pattern for a large distance. Overall, the
use of the measured antenna pattern results in the best fit for the measured
data.

###  V-D Analysis of Shadowing Components from Measurement

(a) 30 m height (b) 90 m height

Figure 11: Shadowing component from measurements and fitted curves to Gaussian
and skewed Gaussian distribution, at (a) 30 m, and (b) 90 m UAV altitude.

After we derive the two-ray path loss model, we can extract the shadowing
component by subtracting the path loss model from measured RSRP using (9), as
shown in Fig. 11. The shadowing component is known to follow a Gaussian
distribution, and the measured shadowing distributions for different UAV
heights are compared to the fitted curves. It is observed that the measured
shadowing distributions can be modeled using a Gaussian distribution, though
there are slight deviations. In particular, the measured distributions exhibit
asymmetry with a heavier left tail compared to the symmetric Gaussian
distribution. To achieve a better fit, an alternative approach is to use a
skewed Gaussian (normal) distribution, which allows for introducing a desired
level of skewness to the distribution [27]. The probability density function
(PDF) of the skewed Gaussian distribution can be expressed as

$\displaystyle
f(x)=2\phi\left(\frac{x-\xi}{\omega}\right)\Phi\left(\alpha\left(\frac{x-\xi}{\omega}\right)\right),$
(38)

where $\phi(\cdot)$, $\Phi(\cdot)$ indicates the PDF and the CDF of Gaussian
distribution, respectively. The parameter $\alpha$ in (38) decides the
skewness of the distribution. If $\alpha$ is a positive real value, it gives
right-skewness, while left-skewness is introduced by a negative real value. In
addition, the mean, the standard deviation of the shadowing, left-skewed
Gaussian parameter $\alpha$, and normalized mean squared error (NMSE) of model
fittings for all heights are listed in Table II. Note that the optimal
$\alpha$ is decided by minimizing NMSE. It shows that the distributions as
well as the value of variances in different heights are similar, and we can
assume a stationary process in spatial data.

TABLE II: Mean and standard deviation of the shadowing component for both measurement and Gaussian fitting curve, and NMSE of models fitting. UAV height | Mean |  Standard deviation | $\alpha$ |  NMSE (Gaussian) |  NMSE (skewed)
---|---|---|---|---|---
$30$ m |  $-4.26$ dB |  $7.14$ dB | -2.13 | 0.0314 | 0.0027
$50$ m |  $-4.57$ dB |  $6.45$ dB | -2.26 | 0.0370 | 0.0020
$70$ m |  $-0.34$ dB |  $6.53$ dB | -2.57 | 0.0437 | 0.0036
$90$ m |  $-0.32$ dB |  $6.90$ dB | -2.08 | 0.0338 | 0.0028
$110$ m |  $-0.27$ dB |  $6.83$ dB | -2.27 | 0.0302 | 0.0022

##  VI Numerical Results on 3D Signal Interpolation

In this section, we will first study the horizontal, vertical, and finally 3D
correlation in the measured data. We will use the 3D correlation to calculate
the semi-variogram, which will subsequently be used to analyze the 3D
interpolation accuracy for various scenarios.

###  VI-A Analysis of Correlation Function from Measurement

####  VI-A1 Horizontal distance correlation

In this subsection, we analyze the spatial correlation using the AERPAW
datasets available at [1]. We obtain correlation functions between two
different 3D locations by using measurements at different heights, and we use
exponential and bi-exponential functions to model the correlations as
discussed earlier. The mean and standard deviation values obtained by
statistical analysis in Section V-D and the measured RSRP values are utilized
in calculating the correlations.

(a) Correlation of horizontal distance.

(b) Correlation of vertical distance.

Figure 12: Correlation of RSRP measurements with horizontal and vertical
distance and curve fits with exponential and bi-expontential models.

We analyze the spatial correlation depending on the horizontal distance
($d_{\rm h}$) with a zero vertical distance ($d_{\rm v}$) by using the
experiment dataset. Since our experiments fix the height of the drone for a
specific flight, the vertical distance between the samples in the same flight
is zero. The analysis of the correlation by the horizontal distance is
performed by following steps:

  1. i.

Calculate the correlation among all samples in a flight, excluding the samples
during the take-off and landing periods.

  2. ii.

Sort the correlation from step (i) according to the horizontal distance
between the sample pairs. This will ensure that the correlations are arranged
in increasing order based on the horizontal distance.

  3. iii.

Average the correlations every 2 m. Start from the smallest horizontal
distance and group the correlations within a 2 m interval. Calculate the
average correlation for each interval. Repeat this process for subsequent 2 m
intervals until covering all the correlations.

  4. iv.

Perform steps (i)-(iii) iteratively for each height (30 m, 50 m, 70 m, 90 m,
110 m). Then, we have correlations for each individual height.

  5. v.

Average the correlation for every distance over all the heights. Take the
correlations obtained in step (iv) for each height and distance, and compute
the average correlation value across all heights for that specific distance.

The correlation between two samples $w_{i}$, $w_{j}$ is calculated by

$\displaystyle
R_{i,j}=\frac{(w_{i}-\nu_{i})(w_{j}-\nu_{j})}{\sigma_{w,i}\sigma_{w,j}},$ (39)

where $\nu$, $\sigma_{w}$ denote the mean and the standard deviation of the
sample, which can be obtained from Table II. The obtained correlation function
and fitted curves are shown in Fig. 12a. It is observed that the correlation
is rapidly decayed as the horizontal distance increases. Although the
correlation is generally modeled by an exponential function (also known as the
Gudmundson model) [10], the bi-exponential model [28] fits better than the
exponential model for our measurements, which is given as

$\displaystyle R(d_{\rm h})=ae^{-b_{1}d_{\rm h}}+(1-a)e^{-b_{2}d_{\rm h}},$
(40)

where $b_{1}$, $b_{2}$ are fitting parameters. We also observe that the
correlation distance is 4.5 m when the correlation is 0.5.

####  VI-A2 Vertical distance correlation

TABLE III: Correlations between different UAV heights. |  $30$ m |  $50$ m |  $70$ m |  $90$ m |  $110$ m
---|---|---|---|---|---
$30$ m | $1.003$ | $0.247$ | $0.080$ | $-0.024$ | $-0.040$
$50$ m | $0.247$ | $1.004$ | $0.214$ | $0.057$ | $-0.002$
$70$ m | $0.080$ | $0.214$ | $1.003$ | $0.307$ | $0.172$
$90$ m | $-0.024$ | $0.057$ | $0.307$ | $1.001$ | $0.409$
$110$ m | $-0.040$ | $-0.002$ | $0.172$ | $0.409$ | $1.012$

We calculate the vertical distance correlation with a zero horizontal distance
from measurements which is opposite to the above subsection. Since the
trajectory of the UAV for flights at different heights is designed to be
identical (see Fig. 4), we can obtain samples of the same 2D location
(latitude, longitude) with different vertical distances. For example, if we
want to obtain 20 m vertical distance samples, we can use the dataset from the
30 m and 50 m UAV flights and pick two samples from any overlapped trajectory
(one from the 30 m height, the other from the 50 m height). The analysis of
the correlation by the vertical distance is conducted by following steps:

  1. i.

Choose two different height measurements datasets, such as the datasets from
the 30 m and 50 m UAV flights;

  2. ii.

Remove data where the two trajectories are not fully overlapped, using a
threshold of $d_{\rm h}>3$ m. This ensures that we have data points with the
same location across the trajectories;

  3. iii.

Calculate the correlations between the two samples with the same location
across the trajectories. Compute the correlation coefficient for each pair of
samples and average them out. This will give you the correlation for a
specific vertical distance (e.g., 20 m) between the two heights.

  4. iv.

Repeat steps (i) to (iii) iteratively for pairs of measurements at different
heights. For example, we can calculate correlations for the 50 m and 70 m
flights, 70 m and 90 m flights, and so on.

In step (ii), we exclude the samples that the trajectory is undesirably not
overlapped by checking GPS readings. The correlations between different pairs
of flights are listed in Table III. We also present the obtained correlation
function from Table III and the fitted curve in Fig. 12b. It is observed that
the correlation function based on the vertical distance fits best with the
exponential model, which is expressed as

$\displaystyle R(d_{\rm v})=e^{-\frac{d_{\rm v}}{d_{\rm cor}}\text{ln}(2)},$
(41)

where the correlation distance is given by $d_{\rm cor}=11.24$ m.

####  VI-A3 3D distance correlation

Figure 13: The 3D distance correlation and fitted curve by the proposed model where vertical distance ($d_{\rm v}$) = 20, 40 m. TABLE IV: Fitting parameters $b_{1}$, $b_{2}$ in (42) depending on the vertical distance. $d_{\rm v}$ | $b_{1}$ | $b_{2}$
---|---|---
0 m | 0.02815 | 0.2474
20 m | 0.05988 | 0.03574
40 m | 0.06998 | 0.045

To analyze the correlation when both horizontal distance and vertical distance
are considered, we can process the dataset obtained from flights at two
different heights. By comparing the measurements from these flights, you can
determine the correlation between two different 3D coordinate locations. The
processing steps for obtaining correlation with 20 m vertical distance are as
follows:

  1. i.

Choose a pair of measurement datasets where the height difference is 20 m. For
example, select the dataset from the 30 m height flight and the dataset from
the 50 m height flight.

  2. ii.

Calculate the correlation between a sample from one height (e.g., 30 m) and a
sample from the other height (e.g., 50 m) across all the samples in the
datasets.

  3. iii.

Sort the correlation from step (ii) by the horizontal distance and average the
correlations for every 2 m of horizontal distance.

  4. iv.

Repeat steps (i) to step (iii) iteratively by different pairs of the
measurement datasets of the height. For example, you can repeat the analysis
with the dataset from the 50 m height flight and the dataset from the 70 m
height flight.

By performing this iterative analysis for different pairs of measurement
datasets with varying height differences, we can obtain the correlation values
that capture the relationship between joint horizontal and vertical distances.
This analysis helps in understanding how the signal strength correlation
varies with changes in both horizontal and vertical distances, providing
insights into the spatial characteristics of the wireless channel.

The 3D distance correlation results with 20 m and 40 m vertical distances are
shown in Fig. 13. We model and fit the correlation of joint horizontal and
vertical distance by combining the correlation functions of the horizontal and
the vertical distance in (40), (41). The proposed correlation model in 3D
space is expressed as

$\displaystyle R(d_{\rm v},d_{\rm h})=e^{-\frac{d_{\rm v}}{d_{\rm
cor}}\text{ln}(2)}\left(ae^{-b_{1}d_{\rm h}}+(1-a)e^{-b_{2}d_{\rm h}}\right),$
(42)

where $a=0.3$, and $b_{1}$, $b_{2}$ are tuning parameters. Note that when
$d_{\rm h}=0$, the model is the same as (40), while when $d_{\rm v}=0$, the
model is equivalent to (41). The fitted values of $b_{1}$, $b_{2}$ depending
on the vertical distance ($d_{\rm v}$) are listed in Table IV.

###  VI-B Analysis of Semi-variogram

Figure 14: The measured semi-variogram and analysis by the correlation
functions in (40) and (41), with respect to horizontal distance (bottom
x-axis), and vertical distance (top x-axis).

In Section III-A, we introduce earlier the concept of semi-variogram in (11)
and derive the relation to the correlation function in (III-A). We analyze the
semi-variogram by measurements results in Fig. 14 with respect to both the
horizontal distance and vertical distance. The measurement results are
directly obtained by the definition of the semi-variogram in (11) and the
analysis results come from the correlation function in (42) which is then used
in (III-A). The measurements and our analysis from (III-A) are closely
overlapped for both distance conditions.

###  VI-C Performance Evaluation with Kriging

In this subsection, we evaluate the 3D interpolation performance of the
Kriging technique described in Section III-B using the measurement dataset. We
adopt cross-validation-based root mean square error (RMSE) evaluation [14],
which compares the predicted RSRP with the measured RSRP to observe the error.
In particular, the RMSE for performance evaluation can be expressed as

$\displaystyle\mathsf{RMSE}=\sqrt{\frac{1}{N_{0}}\sum^{N_{0}}_{i}\left(\hat{r}(l^{\rm
uav}_{0,i})-r(l^{\rm uav}_{0,i})\right)^{2}},$ (43)

where $N_{0}$ denotes the number of samples for prediction. In our evaluation,
the 30 m height measurement samples are predicted by 30 m, 50 m, and 70 m
height measurement datasets. The cross-validation-based evaluation is
conducted by following steps:

  1. i.

Randomly select $M$ samples from the measurement dataset to use for the
prediction. These samples will serve as the training set.

  2. ii.

Randomly select $N_{0}$ samples from the 30 m measurement dataset as the
validation set for cross-validation.

  3. iii.

Use the Kriging technique described in Section III-B to predict the RSRP
values for the $N_{0}$ validation samples based on the $M$ training samples.

  4. iv.

Calculate RMSE between the predicted RSRP values and the actual measured RSRP
values for the $N_{0}$ validation samples. The RMSE is calculated using (43).

  5. v.

Repeat steps (i) to (iv) iteratively for a large number of times, such as
10,000 iterations and calculate the median for the RMSE values obtained from
the iterations. The median value represents the overall prediction performance
of the Kriging technique.

Figure 15: The snapshot of the randomly chosen samples for Kriging.

(a) Prediction by 30 m height measurement.

(b) Prediction by 50 m height measurement.

(c) Prediction by 70 m height measurement.

Figure 16: The RMSE of the prediction at 30 m UAV altitude from measurements
at 30 m, 50 m, and 70 m UAV heights by Kriging, which is benchmarked by the
perfect path loss estimation.

In step (ii), after randomly selecting $M$ samples for prediction, exclude
those samples from the dataset chosen for cross-validation. This ensures that
the samples used for prediction are not used for validation. In addition, when
we predict a sample by Kriging in step (iii), when predicting a sample using
Kriging, consider only the nearby samples within a certain distance threshold
($r_{0}$). Limit the selection of neighboring samples to those within the
$r_{0}$ radius circle around the target sample. These nearby samples will be
used to predict the RSRP value for the target sample. The snapshot of the
randomly chosen $M$ samples from 50 m height measurement and $N_{0}$ from 30 m
height measurement is described in Fig. 15. The figure depicts the radius
circle $r_{0}$ within which nearby samples are used to predict the target
sample. To provide a benchmark for comparison, we consider the perfect path
loss-based 3D interpolation. In particular, we assume that the BS has perfect
knowledge of the exact path loss and transmit power for all locations. This
represents the ideal condition for prediction without utilizing spatial
correlation. The RMSE by the perfect path loss estimation is equivalent to the
standard deviation of the shadowing component from (9) and (43) as follows:

$\displaystyle\mathsf{RMSE}_{\rm
ple}=\sqrt{\mathbb{E}\left[\left(\hat{r}-r\right)^{2}\right]}=\sqrt{\mathbb{E}\left[w^{2}\right]}=\sigma_{w}.$
(44)

In Fig. 16, the RMSE performance of Kriging using measurements at different
UAV altitudes is presented. The results show that the performance of Kriging
varies depending on the altitude of the measurements used for prediction. When
utilizing the 30 m and 50 m height measurement data for prediction, Kriging
outperforms the perfect path loss estimation. This indicates that Kriging can
leverage the spatial correlation present in the highly corrected data to
achieve better prediction accuracy. However, in the case of 70 m height
measurement, the perfect path loss estimation performs better than Kriging.
This suggests that the correlation at a vertical distance of 60 m is too low
to accurately predict using Kriging.

In Fig. 16a, it is observed that the RMSE generally decreases as the number of
samples used for prediction ($N$) increases. However, when $N$ exceeds 250,
the performance of Kriging with an $r_{0}$ value of 200 m is the worst among
the three different r0 values considered. This indicates that while a larger
number of samples can improve performance, adding low-correlated samples can
degrade the prediction accuracy. It is important to strike a balance and
choose an appropriate number of samples ($M$) and radius ($r_{0}$).

Furthermore, in Fig. 16b, the RMSE initially decreases and then increases for
$r_{0}$ values of 70 m, 100 m, and 200 m. This suggests that if the
correlation between samples is not sufficiently high, increasing the number of
samples may not necessarily lead to improved performance. It highlights the
importance of considering both the number of samples and the correlation when
determining the optimal parameters for Kriging prediction.

In conclusion, the choice of the number of samples ($M$) and radius ($r_{0}$)
is crucial for achieving accurate predictions using Kriging. Utilizing a
larger number of highly correlated samples can improve performance, while
including low-correlated samples or selecting an inappropriate radius can
degrade the prediction accuracy.

###  VI-D 3D Interpolation by Kriging

(a) Top view. (b) 3D view.

Figure 17: The 3D radio map from Kriging interpolation using the measurements
at 30 m and 50 m UAV altitudes.

Fig. 17 displays the generated 3D radio map of RSRP using the Kriging
interpolation technique with the available measurement data at 30 m and 50 m
heights. The map provides a visual representation of the RSRP distribution in
the 3D space. The dome shape of the 3D radio map provides valuable insights
into monitoring the signal leakage in the three-dimensional volume of the RDZ.
By examining the map, one can observe the spatial variations and signal
strength levels within the monitored area. The dense 3D radio map obtained
through Kriging interpolation enables efficient analysis and decision-making
related to signal monitoring, interference management, and overall RF planning
within the monitored area. In particular, a spectrum monitoring engine (SME)
can estimate the received signal strength from each signal served within the
RDZ on the surface of the dome. Subsequently, interference to sensitive
receivers outside of the RDZ can be extrapolated, and if exceed a threshold,
interfering signal services in the RDZ can take action (e.g. rescheduling to a
different band or reducing power).

##  VII Conclusion

In this paper, we introduce the RDZ concept which efficiently manages and
controls the spectrum usage by monitoring the signal occupancy and leakage in
a real-time fashion. To monitor the signal leakage from an area, we need to
develop a radio map of signal power surrounding the area, which is more
challenging when considering a 3D space. We propose a signal power
interpolation method in the 3D volume that uses Kriging. The correlation model
between two different 3D locations is designed and the semi-variogram is
defined and analyzed. In addition, we study the proposed 3D Kriging
interpolation using an experimental dataset provided by the NSF AERPAW
platform. We fit path loss and shadowing models to the RSRP measurements and
study the performance of the Kriging interpolation technique for various
scenarios. Our results show that significant gains are possible in received
power estimation accuracy by utilizing the 3D correlation of the data when
compared with using only a path loss based power estimation.

## References

  * [1] S. J. Maeng, O. Ozdemir, I. Guvenc, M. Sichitiu, and R. Dutta, “LTE I/Q Measurement by AERPAW Platform for Air-to-Ground Propagation Modeling,” IEEE Dataport, 2022. [Online]. Available: https://dx.doi.org/10.21227/0p43-0d72
  * [2] M. Zheleva, C. R. Anderson, M. Aksoy, J. T. Johnson, H. Affinnih, and C. G. DePree, “Radio Dynamic Zones: Motivations, challenges, and opportunities to catalyze spectrum coexistence,” _IEEE Commun. Mag._ , 2023.
  * [3] S. J. Maeng, I. Güvenç, M. Sichitiu, B. Floyd, R. Dutta, T. Zajkowski, O. Ozdemir, and M. Mushi, “National radio dynamic zone concept with autonomous aerial and ground spectrum sensors,” in _IEEE Int. Conf. Commun. Workshops (ICC Workshops)_ , Seoul, Korea, Republic of, May 2022, pp. 687–692.
  * [4] S. Tschimben, A. Aradhya, G. Weihe, M. Lofquist, A. Pollak, W. Farah, D. DeBoer, and K. Gifford, “Testbed for Radio Astronomy Interference Characterization and Spectrum Sharing Research,” in _IEEE Aerosp. Conf._ , Big Sky, MT, USA, Mar. 2023, pp. 1–16.
  * [5] “Spectrum innovation initiative: National radio dynamic zones (sii-nrdz),” NSF Program Solicitation, June 2022. [Online]. Available: https://www.nsf.gov/pubs/2022/nsf22579/nsf22579.htm
  * [6] H. Zou, M. Jin, H. Jiang, L. Xie, and C. J. Spanos, “Winips: Wifi-based non-intrusive indoor positioning system with online radio map construction and adaptation,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 12, pp. 8118–8130, Dec. 2017.
  * [7] 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.
  * [8] F. Graziosi and F. Santucci, “A general correlation model for shadow fading in mobile radio systems,” _IEEE Commun. Lett._ , vol. 6, no. 3, pp. 102–104, Mar. 2002.
  * [9] S. S. Szyszkowicz, H. Yanikomeroglu, and J. S. Thompson, “On the feasibility of wireless shadowing correlation models,” _IEEE Trans. Veh. Technol._ , vol. 59, no. 9, pp. 4222–4236, Nov. 2010.
  * [10] M. Gudmundson, “Correlation model for shadow fading in mobile radio systems,” _Electronics letters_ , vol. 23, no. 27, pp. 2145–2146, 1991.
  * [11] R. He, Z. Zhong, B. Ai, and C. Oestges, “Shadow fading correlation in high-speed railway environments,” _IEEE Trans. Veh. Technol._ , vol. 64, no. 7, pp. 2762–2772, Jul. 2015.
  * [12] H. Braham, S. B. Jemaa, G. Fort, E. Moulines, and B. Sayrac, “Fixed rank kriging for cellular coverage analysis,” _IEEE Trans. Veh. Technol._ , vol. 66, no. 5, pp. 4212–4222, May 2017.
  * [13] K. Sato and T. Fujii, “Kriging-based interference power constraint: Integrated design of the radio environment map and transmission power,” _IEEE Trans. Cogn. Commun. Netw._ , vol. 3, no. 1, pp. 13–25, Mar. 2017.
  * [14] K. Sato, K. Suto, K. Inage, K. Adachi, and T. Fujii, “Space-frequency-interpolated radio map,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 1, pp. 714–725, Jan. 2021.
  * [15] A. Al-Hourani and K. Gomez, “Modeling Cellular-to-UAV Path-Loss for Suburban Environments,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 1, pp. 82–85, Feb. 2018.
  * [16] J. Holis and P. Pechac, “Elevation dependent shadowing model for mobile communications via high altitude platforms in built-up areas,” _IEEE Trans. Antennas Propag._ , vol. 56, no. 4, pp. 1078–1084, Apr. 2008.
  * [17] M. Simunek, F. P. Fontán, and P. Pechac, “The UAV Low Elevation Propagation Channel in Urban Areas: Statistical Analysis and Time-Series Generator,” _IEEE Trans. Antennas Propag._ , vol. 61, no. 7, pp. 3850–3858, Jul. 2013\.
  * [18] S. J. Maeng, İ. Güvenç, M. L. Sichitiu, and O. Ozdemir, “Out-of-zone signal leakage sensing in radio dynamic zones,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , Seoul, Korea, May 2022.
  * [19] T. Kidd, “National radio quiet and dynamic zones,” CHIPS – The Department of Navy’s Information Technology Magazine, Apr.-June 2018. [Online]. Available: https://www.doncio.navy.mil/CHIPS/ArticleDetails.aspx?ID=10299
  * [20] N. R. Chopde and M. Nichat, “Landmark based shortest path detection by using A* and Haversine formula,” _International Journal of Innovative Research in Computer and Communication Engineering_ , vol. 1, no. 2, pp. 298–302, Apr. 2013.
  * [21] W. C. Jakes and D. C. Cox, _Microwave mobile communications_. Wiley-IEEE press, 1994.
  * [22] N. Cressie, _Statistics for spatial data_. John Wiley & Sons, 2015.
  * [23] Octane Wireless, “SA-1400-5900 Data Sheet.” [Online]. Available: https://www.octanewireless.com/product/sa-1400-5900-tri-band-stub-antenna/
  * [24] Mobile Mark, Inc, “RM-WB1 Series Radiation Pattern.” [Online]. Available: https://www.mobilemark.com/product/rm-wb1/
  * [25] S. J. Maeng, O. Ozdemir, İ. Güvenç, M. Sichitiu, R. Dutta, and M. Mushi, “AERIQ: SDR-Based LTE I/Q Measurement and Analysis Framework for Air-to-Ground Propagation Modeling,” in _IEEE Aerosp. Conf._ , Big Sky, MT, USA, Mar. 2023, pp. 1–11.
  * [26] S. J. Maeng, M. A. Deshmukh, I. Güvenç, A. Bhuyan, and H. Dai, “Interference analysis and mitigation for aerial IoT considering 3D antenna patterns,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 1, pp. 490–503, Jan. 2021.
  * [27] C. K. Sung, S. Li, M. Hedley, N. Nikolic, and W. Ni, “Skew log-normal channel model for indoor cooperative localization,” in _Proc. IEEE Int. Symp. Pers., Indoor, Mobile Radio Commun. (PIMRC)_ , Montreal, Canada, Oct. 2017, pp. 1–5.
  * [28] L. Liu, C. Tao, D. W. Matolak, T. Zhou, and H. Chen, “Investigation of shadowing effects in typical propagation scenarios for high-speed railway at 2350 mhz,” _International Journal of Antennas and Propagation_ , vol. 2016, Oct. 2016.