5G-NR

    5G New Radio
3GPP

    3rd Generation Partnership Project
ABS

    aerial base station
AC

    address coding
ACF

    autocorrelation function
ACR

    autocorrelation receiver
ADC

    analog-to-digital converter
AIC

    Analog-to-Information Converter
AIC

    Akaike information criterion
ARIC

    asymmetric restricted isometry constant
ARIP

    asymmetric restricted isometry property
ARQ

    Automatic Repeat Request
AUB

    asymptotic union bound
AWGN

    Additive White Gaussian Noise
AWGN

    additive white Gaussian noise
PSK

    asymmetric PSK
AWRICs

    asymmetric weak restricted isometry constants
AWRIP

    asymmetric weak restricted isometry property
BCH

    Bose, Chaudhuri, and Hocquenghem
BCHSC

    BCH based source coding
BEP

    bit error probability
BFC

    block fading channel
BG

    Bernoulli-Gaussian
BGG

    Bernoulli-Generalized Gaussian
BPAM

    binary pulse amplitude modulation
BPDN

    Basis Pursuit Denoising
BPPM

    binary pulse position modulation
BPSK

    Binary Phase Shift Keying
BPZF

    bandpass zonal filter
BSC

    binary symmetric channels
BU

    Bernoulli-uniform
BER

    bit error rate
BS

    base station
BW

    BandWidth
BLLL

     binary log-linear learning
CP

    Cyclic Prefix
CDF

    cumulative distribution function
CDF

    Cumulative Distribution Function
CDF

    cumulative distribution function
CCDF

    complementary cumulative distribution function
CCDF

    complementary CDF
CCDF

    complementary cumulative distribution function
CD

    cooperative diversity
CDMA

    Code Division Multiple Access
ch.f.

    characteristic function
CIR

    channel impulse response
CoSaMP

    compressive sampling matching pursuit
CR

    cognitive radio
CS

    compressed sensing
CS

    Compressed sensing
CS

    compressed sensing
CSI

    channel state information
CCSDS

    consultative committee for space data systems
CC

    convolutional coding
COVID-19

    Coronavirus disease
DAA

    detect and avoid
DAB

    digital audio broadcasting
DCT

    discrete cosine transform
DFT

    discrete Fourier transform
DR

    distortion-rate
DS

    direct sequence
DS-SS

    direct-sequence spread-spectrum
DTR

    differential transmitted-reference
DVB-H

    digital video broadcasting – handheld
DVB-T

    digital video broadcasting – terrestrial
DL

    DownLink
DSSS

    Direct Sequence Spread Spectrum
DFT-s-OFDM

    Discrete Fourier Transform-spread-Orthogonal Frequency Division Multiplexing
DAS

    Distributed Antenna System
DNA

    DeoxyriboNucleic Acid
EC

    European Commission
EED

    exact eigenvalues distribution
EIRP

    Equivalent Isotropically Radiated Power
ELP

    equivalent low-pass
eMBB

    Enhanced Mobile Broadband
EMF

    ElectroMagnetic Field
EU

    European union
EI

    Exposure Index
eICIC

    enhanced Inter-Cell Interference Coordination
FC

    fusion center
FCC

    Federal Communications Commission
FEC

    forward error correction
FFT

    fast Fourier transform
FH

    frequency-hopping
FH-SS

    frequency-hopping spread-spectrum
FS

    Frame synchronization
FS

    frame synchronization
FDMA

    Frequency Division Multiple Access
GA

    Gaussian approximation
GF

    Galois field
GG

    Generalized-Gaussian
GIC

    generalized information criterion
GLRT

    generalized likelihood ratio test
GPS

    Global Positioning System
GMSK

    Gaussian Minimum Shift Keying
GSMA

    Global System for Mobile communications Association
GS

    ground station
GMG

     Grid-connected MicroGeneration
HAP

    high altitude platform
HetNet

    Heterogeneous network
IDR

    information distortion-rate
IFFT

    inverse fast Fourier transform
IHT

    iterative hard thresholding
i.i.d.

    independent, identically distributed
IoT

    Internet of Things
IR

    impulse radio
LRIC

    lower restricted isometry constant
LRICt

    lower restricted isometry constant threshold
ISI

    intersymbol interference
ITU

    International Telecommunication Union
ICNIRP

    International Commission on Non-Ionizing Radiation Protection
IEEE

    Institute of Electrical and Electronics Engineers
ICES

    IEEE international committee on electromagnetic safety
IEC

    International Electrotechnical Commission
IARC

    International Agency on Research on Cancer
IS-95

    Interim Standard 95
KPI

    Key Performance Indicator
LEO

    low earth orbit
LF

    likelihood function
LLF

    log-likelihood function
LLR

    log-likelihood ratio
LLRT

    log-likelihood ratio test
LoS

    Line-of-Sight
LRT

    likelihood ratio test
LWRIC

    lower weak restricted isometry constant
LWRICt

    LWRIC threshold
LPWAN

    Low Power Wide Area Network
LoRaWAN

    Low power long Range Wide Area Network
NLoS

    Non-Line-of-Sight
Li-Fi

    light-fidelity
LED

    light emitting diode
LABS

    LoS transmission with each ABS
NLABS

    NLoS transmission with each ABS
MB

    multiband
MC

    macro cell
MDS

    mixed distributed source
MF

    matched filter
m.g.f.

    moment generating function
MI

    mutual information
MIMO

    Multiple-Input Multiple-Output
MISO

    multiple-input single-output
MJSO

    maximum joint support cardinality
ML

    maximum likelihood
MMSE

    minimum mean-square error
MMV

    multiple measurement vectors
MOS

    model order selection
${M}$-PSK

     $M$-ary phase shift keying
${M}$-PSK

     $M$-ary asymmetric PSK
MP

     multi-period
MINLP

    mixed integer non-linear programming
$M$-QAM

     $M$-ary quadrature amplitude modulation
MRC

    maximal ratio combiner
MSO

    maximum sparsity order
M2M

    Machine-to-Machine
MUI

    multi-user interference
mMTC

    massive Machine Type Communications
mm-Wave

    millimeter-wave
MP

    mobile phone
MPE

    maximum permissible exposure
MAC

    media access control
NB

    narrowband
NBI

    narrowband interference
NLA

    nonlinear sparse approximation
NLOS

    Non-Line of Sight
NTIA

    National Telecommunications and Information Administration
NTP

    National Toxicology Program
NHS

    National Health Service
LOS

    Line of Sight
OC

    optimum combining
OC

    optimum combining
ODE

    operational distortion-energy
ODR

    operational distortion-rate
OFDM

    Orthogonal Frequency-Division Multiplexing
OMP

    orthogonal matching pursuit
OSMP

    orthogonal subspace matching pursuit
OQAM

    offset quadrature amplitude modulation
OQPSK

    offset QPSK
OFDMA

    Orthogonal Frequency-division Multiple Access
OPEX

    Operating Expenditures
OQPSK/PM

    OQPSK with phase modulation
PAM

    pulse amplitude modulation
PAR

    peak-to-average ratio
PDF

    probability density function
PDF

    probability density function
PDF

    probability distribution function
PDP

    power dispersion profile
PMF

    probability mass function
PMF

    probability mass function
PN

    pseudo-noise
PPM

    pulse position modulation
PRake

    Partial Rake
PSD

    power spectral density
PSEP

    pairwise synchronization error probability
PSK

    phase shift keying
PD

    power density
$8$-PSK

     $8$-phase shift keying
PPP

    Poisson point process
PCP

    Poisson cluster process
FSK

    Frequency Shift Keying
QAM

    Quadrature Amplitude Modulation
QPSK

    Quadrature Phase Shift Keying
OQPSK/PM

    OQPSK with phase modulator
RD

    raw data
RDL

    ”random data limit”
RIC

    restricted isometry constant
RICt

    restricted isometry constant threshold
RIP

    restricted isometry property
ROC

    receiver operating characteristic
RQ

    Raleigh quotient
RS

    Reed-Solomon
RSSC

    RS based source coding
r.v.

    random variable
R.V.

    random vector
RMS

    root mean square
RFR

    radiofrequency radiation
RIS

    Reconfigurable Intelligent Surface
RNA

    RiboNucleic Acid
RRM

    Radio Resource Management
RUE

    reference user equipments
RAT

    radio access technology
RB

    resource block
SA-Music

    subspace-augmented MUSIC with OSMP
SC

    small cell
SCBSES

    Source Compression Based Syndrome Encoding Scheme
SCM

    sample covariance matrix
SEP

    symbol error probability
SG

    sparse-land Gaussian model
SIMO

    single-input multiple-output
SINR

    signal-to-interference plus noise ratio
SIR

    signal-to-interference ratio
SISO

    Single-Input Single-Output
SMV

    single measurement vector
SNR

    signal-to-noise ratio
SP

    subspace pursuit
SS

    spread spectrum
SW

    sync word
SAR

    specific absorption rate
SSB

    synchronization signal block
SR

    shrink and realign
tUAV

    tethered Unmanned Aerial Vehicle
TBS

    terrestrial base station
uUAV

    untethered Unmanned Aerial Vehicle
PDF

    probability density functions
PL

    path-loss
TH

    time-hopping
ToA

    time-of-arrival
TR

    transmitted-reference
TW

    Tracy-Widom
TWDT

    TW Distribution Tail
TCM

    trellis coded modulation
TDD

    Time-Division Duplexing
TDMA

    Time Division Multiple Access
Tx

    average transmit
UAV

    Unmanned Aerial Vehicle
URIC

    upper restricted isometry constant
URICt

    upper restricted isometry constant threshold
UWB

    ultrawide band
UWB

    Ultrawide band
URLLC

    Ultra Reliable Low Latency Communications
UWRIC

    upper weak restricted isometry constant
UWRICt

    UWRIC threshold
UE

    User Equipment
UL

    UpLink
WiM

    weigh-in-motion
WLAN

    wireless local area network
WM

    Wishart matrix
WMAN

    wireless metropolitan area network
WPAN

    wireless personal area network
WRIC

    weak restricted isometry constant
WRICt

    weak restricted isometry constant thresholds
WRIP

    weak restricted isometry property
WSN

    wireless sensor network
WSS

    Wide-Sense Stationary
WHO

    World Health Organization
Wi-Fi

    Wireless Fidelity
SpaSoSEnc

    sparse source syndrome encoding
VLC

    Visible Light Communication
VPN

    Virtual Private Network
RF

    Radio Frequency
FSO

    Free Space Optics
IoST

    Internet of Space Things
GSM

    Global System for Mobile Communications
2G

    Second-generation cellular network
3G

    Third-generation cellular network
4G

    Fourth-generation cellular network
5G

    Fifth-generation cellular network
gNB

    next-generation Node-B Base Station
NR

    New Radio
UMTS

    Universal Mobile Telecommunications Service
LTE

    Long Term Evolution
QoS

    Quality of Service

#  Resident Population Density-Inspired Deployment of K-tier Aerial Cellular
Network

Ruibo Wang, Mustafa A. Kishk, Member, IEEE and Mohamed-Slim Alouini, Fellow,
IEEE Ruibo Wang and Mohamed-Slim Alouini are with King Abdullah University of
Science and Technology (KAUST), CEMSE division, Thuwal 23955-6900, Saudi
Arabia. Mustafa A. Kishk is with the Department of Electronic Engineering,
National University of Ireland, Maynooth, W23 F2H6, Ireland. (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>slim.alouini@kaust.edu.sa).

###### Abstract

Using Unmanned Aerial Vehicles to enhance network coverage has proven a
variety of benefits compared to terrestrial counterparts. One of the commonly
used mathematical tools to model the locations of the UAVs is stochastic
geometry (SG). However, in the existing studies, both users and UAVs are often
modeled as homogeneous point processes. In this paper, we consider an
inhomogeneous Poisson point process (PPP)-based model for the locations of the
users that captures the degradation in the density of active users as we move
away from the town center. In addition, we propose the deployment of aerial
vehicles following the same inhomogeneity of the users to maximize the
performance. In addition, a multi-tier network model is also considered to
make better use of the rich space resources. Then, the analytical expressions
of the coverage probability for a typical user and the total coverage
probability are derived. Finally, we optimize the coverage probability with
limitations of the total number of UAVs and the minimum local coverage
probability. Finally we give the optimal UAV distribution parameters when the
maximum overall coverage probability is reached.

###### Index Terms:

Coverage probability, urban model, multi-tier UAV network, stochastic
geometry, inhomogeneous Poisson point process.

##  I Introduction

In the next generation mobile network (5G, beyond 5G), UAVs have many
application scenarios [1, 2, 3], among which UAV-aided ubiquitous coverage
becomes an important topic [4]. Because UAVs are easy to deploy, highly
mobile, and have 3D deployment, they are often used to build temporary or
dynamic networks and provide ubiquitous coverage. Especially, UAV is widely
used to relieve the pressure of large crowds gathering in small areas [4, 5].
They are proven to provide reliable system coverage in hot spots and provide
additional system performance [6]. Due to the demand for high rate signals,
multi-tier vertical heterogeneous networks (VHetNet) is proposed to make use
of space resources in city centers [7, 8].

One of the main unanswered questions in the realm of UAV-enabled wireless
networks is where and how high the UAVs should be deployed [9]. The common
assumption in SG-based literature is that the user’s spatial distribution is
homogeneous. Consequently, existing literature typically assumes that the
density of the UAVs is spatially invariant. However, according to recent
studies on resident population densities, a more proper assumption would be
for the density of the users (and consequently the UAVs) drops as their
distance from the town center increases [10]. Analyzing the influence of such
a setup on the wireless network’s performance is the main objective of this
paper. More details on the contributions of this paper are provided later in
Sec. I-B.

###  I-A Related Work

SG is a powerful mathematical method of analyzing communication networks with
irregular topology [11]. Furthermore, the SG framework is suitable for
modeling and analyzing devices in motion, such as UAVs, cars [12], and LEO
satellites [13, 14]. The SG-based analytical results of the network coverage
probability can provide accurate approximations to the actual network [15,
16]. Next, the authors in [17] proposed an air-to-ground line-of-sight (LoS)
probability model suitable for town centers. In this model, the probability of
the UAV being blocked by the building decreases with the increase of the
elevation angle of the UAV to the typical user. This model divides UAVs into
LoS UAVs and non-line-of-sight (NLoS) UAVs. Since the model is related to
density, area, and height of building [18], it is suitable for various
scenarios. Based on the LoS probability model, there has been some literature
on UAV networking in town centers [19, 6].

In the existing research, some resident population density models have been
considered. Different user distributions in several urban environments are
proposed in [20]. A disjoint clustered model for large resident population
density is set up in [21]. A central model is provided in [10] and is adopted
in this paper. In the central model, user density decreases with the distance
from the user to the center. However, the above articles pay more attention to
the modeling of users, while the UAVs are simply deployed. UAVs are deployed
as a homogeneous PPP in [10, 21], while the locations of UAVs are determined
by clustering in [20]. Therefore, the deployment of UAVs is also worth
exploring. However, with regard to analyzing downlink network coverage
performance, changing the distribution of UAVs brings much more difficulty in
technical derivation than changing the distribution of users. Given that UAVs
form a homogeneous PPP, the downlink coverage performance of users at any
location is the same. Nevertheless, when the density of the UAV is not
constant, the distributions of the distance between the serving UAV and the
interfering UAV to the user are different for the users at different
locations, which makes the analysis challenging.

To effectively utilize the deployable space of UAVs, developing the vertical
deployment mode of UAVs is also worth studying, in addition to designing the
horizontal distribution of UAVs. Based on the SG framework, authors in [8, 9,
22, 23, 24, 25] have put forward multi-tier VHetNets consisting of ground base
stations, UAVs and high altitude platforms and low earth orbit
(LEO)-satellites. The above researches all introduced the concept of
association probability to describe the probability of users choosing a
communication device in a specific tier (instead of other tiers) to provide
services. Unfortunately, the above analytical framework is unsuitable for our
study because the UAVs are not uniformly distributed in our paper. Designing a
different method to obtain the association probability is another challenge.

###  I-B Contribution

The contributions of this paper can be summarized as follows:

  * •

We study a resident population density-inspired model of the urban area. The
density of users decreases with the distance to the town center. UAVs follow a
similar distribution to the distribution of users and are deployed at
different altitudes with different densities.

  * •

We derive the analytical result of coverage probability under the specific
model and prove that it is consistent with the Monte-Carlo simulation. In
addition, the existing coverage probability analysis framework is extended to
data rate and energy efficiency.

  * •

The coverage performance of multi-tier networks and single-tier networks are
compared. We also compare the coverage performance of the population density-
inspired distribution and the homogeneous distribution of UAV.

  * •

By adjusting the distribution of UAVs in each tier, we optimize the coverage
probability under different user distributions. Furthermore, remarks on the
parameter design criterion for UAV distribution are given.

##  II System Model

###  II-A Network Model

Figure 1: Illustration of the system model.

As shown in Fig. 1, we consider a scenario in which ground users are
distributed according to an inhomogeneous PPP, which is inspired by the
resident population distribution model proposed in [10]. We assume that the
UAVs in the VHetNet are deployed based on the ground users density and, hence,
their locations also follow an inhomogeneous PPP. Assuming that the center of
the town is located at the origin, the densities of the users and the UAVs
near the origin are relatively high, while the density goes down as we move
away from the origin. Assume the users are located at the ground, with
horizontal distance $z_{u}$ to the origin, the density distribution of the
users ${\Lambda_{u}}\left(z_{u}\right)$ can be represented as follows,

${\Lambda_{u}}\left(z_{u}\right)={\lambda_{u}}{e^{-{\beta_{u}}z_{u}}},$ (1)

where $\lambda_{u}$ determines the total density of the plane, $\beta_{u}$ is
a measure of homogeneity. When the value of $\beta_{u}$ is large, the users
are spatially condensed at the origin. When $\beta_{u}=0$ the process
degenerates to a homogeneous PPP. Without loss of generality, we focus on a
typical user located on the positive $X$-axis.

We assume that $K$ tiers of UAVs are distributed at a set of some fixed
heights $h_{k}$ independently. Their location distribution of each tier form a
2D inhomogeneous PPP, denoted by
${\Phi_{k}}\mathop{=}\limits^{\Delta}\left\\{{{x_{i,k}}}\right\\}$ where
${{x_{i,k}}}$ refers to the 3D location of UAV $i$ in tier $k$. We prefer
polar coordinates $\left(z_{i},\theta_{i},h_{k}\right)$ to represent
${x_{i,k}}$, where $z_{i}$ is the horizontal distance between the UAV and the
origin, $\theta_{i}$ is the angle between the $X$-axis and the line which
connects the projection of the UAV and the origin. Because designers tend to
place more UAVs in densely populated areas, it is reasonable to assume that
they will be in the same distribution as the users. Thus, in tier $k$, the
density distribution ${\Lambda_{\rm{UAV},k}}$ can be described as,

${\Lambda_{\rm{UAV},k}}\left(z_{i}\right)={\lambda_{k}}{e^{-{\beta_{k}}z_{i}}},$
(2)

where $\lambda_{k}$ and $\beta_{k}$ are parameters of the UAVs in tier k,
which have the same meaning as the users’ parameters in the density
distribution. Furthermore, we assume that each tier of UAVs have the same
transmitting power ${\rho_{k}}$. The quad
${T_{\rm{k}}}=\left\\{{{h_{k}},{\lambda_{k}},{\beta_{k}},{\rho_{k}}}\right\\},k=1,2,...,K$
is used to represent the parameters of the ${k^{th}}$ tier.

###  II-B Channel Model

To model the air-to-ground channel between a user and a UAV, we need to take
into consideration the LoS and NLoS scenarios [17]. Considering a UAV in tier
$k$, given the horizontal distance $z$ between the UAV’s projection on the
ground and the user, the probability of setting up an LoS link between the
typical user and the UAV is [17], [26],

$P_{k}^{{\rm{LoS}}}\left({\rm{z}}\right)=\frac{1}{{1+a\exp\left({-b\left({\frac{{180}}{\pi}{{\tan}^{-1}}\left({\frac{{{h_{k}}}}{z}}\right){-a}}\right)}\right)}},$
(3)

where $a$ and $b$ are environment-dependent parameters. From the perspective
of the typical user, the inhomogeous PPP process corresponding to the $K$-tier
UAVs can be split into two disjoint PPPs, that is
${\Phi_{k}}{=}{\Phi_{{\rm{LoS}},k}}\cup{\Phi_{{\rm{NLoS}},k}}$ and
${\Phi_{{\rm{LoS}},k}}\cap{\Phi_{{\rm{NLoS}},k}}=\phi$, where
${\Phi_{{\rm{LoS}},k}}$ and ${\Phi_{{\rm{NLoS}},k}}$ denote the set of UAVs
which establish LoS and NLoS conditions for the typical user respectively,
$\phi$ is the empty set.

In this article, UAVs with an LoS link to the typical user are abbreviated as
LoS UAVs, while the rest are abbreviated as NLoS UAVs. Therefore, the 3D
VHetNet is split into $2K$ disjoint two-dimensional PPPs, with the density
$P_{k}^{Q}\,{\Lambda_{\rm{UAV},k}}$, where
$Q=\left\\{{{\rm{LoS}},{\rm{NLoS}}}\right\\}$. After LoS and NLoS states are
defined, the channel fading model can be established, which is described by
small-scale fading and large-scale fading.

For small scale fading, we denote channel fading power gains in terms of
independent random variables $G_{\rm{LoS}}$ and $G_{\rm{NLoS}}$, under LoS and
NLoS conditions for the typical user, respectively. In order to represent
several fading scenarios, Nakagami-m fading is experienced with shape
parameters and scale parameters $(m_{\rm{LoS}},\frac{1}{m_{\rm{LoS}}})$ and
$(m_{\rm{NLoS}},\frac{1}{m_{\rm{NLoS}}})$ for LoS and NLoS links,
respectively. As a result, the probability density functions (PDF) of the
power gains $G_{Q}$ is given by [27]

${f_{G_{Q}}}\left(g\right)=\frac{{{m_{Q}}^{{m_{Q}}}{r^{{m_{Q}}-1}}}}{{\Gamma\left({{m_{Q}}}\right)}}{e^{-{m_{Q}}\,g}},$
(4)

where
$\Gamma\left({{m_{Q}}}\right)=\int_{0}^{\infty}{{x^{{m_{Q}}-1}}{e^{-x}}dx}$ is
the Gamma function, $Q=\left\\{\rm{LoS},\rm{NLoS}\right\\}$.

For large scale fading, $\eta_{\rm{LoS}}$ and $\eta_{\rm{NLoS}}$ are mean
additional gain for LoS and NLoS transmissions [17], with
${\eta_{\rm{LoS}}}>{\eta_{\rm{NLoS}}}$ satisfied. Combining small scale and
large scale fading, the received power of the typical user, transmitted by a
UAV in tier k, is given by,

$\displaystyle{S_{k}}\left(r\right)=\left\\{\begin{array}[]{ll}{{\eta_{\rm{LoS}}}{\rho_{k}}{G_{\rm{LoS}}}{r^{-{\alpha_{\rm{LoS}}}}}}&{{\rm{in\,case\,of\,LoS}}}\\\
{\eta_{\rm{NLoS}}}{\rho_{k}}{G_{\rm{NLoS}}}{r^{-{\alpha_{\rm{NLoS}}}}}&{{\rm{in\,case\,of\,NLoS}}}\end{array}.\right.$
(7)

where ${\rho_{k}}$ is the transmission power of UAVs in tier $k$,
$\alpha_{\rm{LoS}}$ and $\alpha_{\rm{NLoS}}$ are path-loss exponents for LoS
and NLoS transmissions, with ${\alpha_{\rm{LoS}}}<{\alpha_{\rm{NLoS}}}$
satisfied, $r$ is the Euclidean distance between the UAV and the typical user,
which is computed by polar coordinates $\left(z_{i},\theta_{i},h_{k}\right)$
of the UAV and the horizontal distance $z_{u}$ from the user to the origin

$r_{{x_{i},k}}=\sqrt{{{({z_{i}}{\rm{cos}}{\theta_{i}}-{z_{u}})}^{2}}+{{({z_{i}}{\rm{sin}}{\theta_{i}})}^{2}}+h_{k}^{2}}.$
(8)

###  II-C Interference

In each tier, the closest LoS and NLoS UAVs are called tagged UAVs. According
to the strongest average received power association strategy [9], the typical
user will associate with the UAV with the strongest average received power
among these $2K$ tagged UAVs. Furthermore, denote the location of the
associated UAV as $x_{o}$. Note that the typical user may not associate with
the closest UAV, since in different tiers, the transmitted power of the UAVs
is different.

In urban areas where UAVs are densely distributed, it is necessary to consider
interference between UAVs. Considering a worst-case scenario, except for the
associated UAV, other LoS and NLoS UAVs interfere with the typical user, which
we refer to as interfering UAVs. Given that the horizontal distance from the
typical user to the origin is $z_{u}$ and the associated UAV is in tier $j$,
the total interference can be expressed as a function of $r_{x_{o}}$, which is
the Euclidean distance between the associated UAV and the typical user,

$\begin{split}{I_{j}}\left({r_{x_{o}},{z_{u}}}\right)&=\sum\limits_{k=1}^{K}\Big{(}\sum\limits_{{x_{i,k}}\in{\Phi_{{\rm{LoS}},k}}\backslash\left\\{{{x_{o}}}\right\\}}{{\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}r_{{x_{i},k}}^{-{\alpha_{{\rm{LoS}}}}}}+\sum\limits_{{x_{i,k}}\in{\Phi_{{\rm{NLoS}},k}}\backslash\left\\{{{x_{o}}}\right\\}}{{\eta_{{\rm{NLoS}}}}{\rho_{k}}{G_{{\rm{NLoS}}}}r_{{x_{i},k}}^{-{\alpha_{{\rm{NLoS}}}}}}\Big{)}.\end{split}$
(9)

As shown in the above formulation, the value of the total interference is
related to $z_{u}$, $r_{x_{o}}$ and the tier where the associated UAV is
located. $r_{x_{o}}$ and the transmission power of tier $j$ determine the
received power from the associated UAV.

###  II-D Performance Analysis

Assuming the associated UAV is located in tier $j$, the instantaneous signal-
to-interference plus noise ratio (SINR) at the typical user is given by the
following equation,

$\displaystyle{\rm{SINR}}=\left\\{\begin{array}[]{ll}{\frac{{{\eta_{{\rm{LoS}}}}{\rho_{j}}{G_{{\rm{LoS}}}}r_{{x_{o}}}^{{}^{-{\alpha_{{\rm{LoS}}}}}}}}{{{I_{j}}\left({{r_{x_{o}}}\left|{{z_{u}}}\right.}\right)+{\sigma^{2}}}}}&{{x_{o}}\in{\Phi_{{\rm{LoS}},k}}}\\\
{\frac{{{\eta_{{\rm{NLoS}}}}{\rho_{j}}{G_{{\rm{NLoS}}}}r_{{x_{o}}}^{{}^{-{\alpha_{{\rm{NLoS}}}}}}}}{{{I_{j}}\left({{r_{{x_{o}}}}\left|{{z_{u}}}\right.}\right)+{\sigma^{2}}}}}&{{x_{o}}\in{\Phi_{{\rm{NLoS}},k}}}\end{array}.\right.$
(12)

where $\sigma^{2}$ is the additive white Gaussian noise (AWGN) power, and
${I_{j}}\left({r_{x_{o}}\left|{{z_{u}}}\right.}\right)$ is the total
interference power.

The reliability of the service provided can be evaluated by the average
performance of SINR. In consequence, the coverage probability, which
represents the probability that the system can provide reliable connections is
defined as the probability that the SINR is greater than a predefined
threshold $\gamma$:

${P^{C}}\mathop{=}\limits^{\Delta}\mathbb{P}\left[{\rm{SINR}>\gamma}\right].$
(13)

##  III Problem Formulation

In this section, our objective is to obtain the analytical expression of
coverage probability. From the definition of coverage probability in (13), we
know that distinguishing the associated UAV and the interfering UAVs is a
prerequisite for computing the coverage probability. Taking an LoS associated
UAV for example, the following steps are used to obtain the analytic
expression for the coverage probability: (i) derive the PDF of the distance
distribution of tagged UAV in each tier $k$, denoted as
${f_{{\rm{RLoS}},k}}\left({r,z_{u}}\right)$, (ii) calculate the probability of
the tagged UAV in tier $k$ being associated with the typical user, defined as
association probability $P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)$, with
${f_{{\rm{RLoS}},k}}\left({r,z_{u}}\right)P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)$
being the PDF of the distance between the associated UAV in tier $k$ and the
typical user, (iii) for a specific distance $r$ between the typical user and
its associated UAV, denote the probability that the SINR is greater than the
threshold $\gamma$ as the conditional coverage probability
${P^{C}}\left({\gamma,{z_{u}}\left|r\right.}\right)$. The average coverage
probability is obtained by taking the expectation of the conditional coverage
probability ${P^{C}}\left({\gamma,r,{z_{u}}}\right)$ with respect to the PDF
of $r$ in step (ii). Steps (i) - (iii) will be explained in Sec. III-B Sec.
III-C, and Sec. III-D, respectively.

###  III-A Nearest Interfering UAVs

A clear understanding of the location range of interfering UAVs is necessary
when analyzing tagged or associated drone distribution. Obviously, for tier
$k$, the nearest interfering UAV should locate at a distance larger than
$h_{k}$ for the typical user. Another possible lower bound of the distance for
an interfering UAV is to ensure a lower average receiving power than the
associated UAV. In the following lemmas, the distance between the typical user
and the nearest interfering UAVs is given.

######  Lemma 1.

Given that the typical user is associated with a LoS UAV located at distance
$r$ in tier $j$, the closest interfering LoS and NLoS UAV in tier $k$ are at
least at distances ${d_{{\rm{LoS}}{\rm{-}}{\rm{LoS}},j,k}}\left(r\right)$ and
${d_{{\rm{LoS}}{\rm{-}}\rm{NLoS},j,k}}\left(r\right)$ , given by

${d_{{\rm{LoS-
LoS}},j,k}}\left(r\right)=\max\left\\{{{h_{k}},{{\left({\frac{{{\rho_{k}}}}{{{\rho_{j}}}}}\right)}^{\frac{1}{{{\alpha_{{\rm{LoS}}}}}}}}r}\right\\},$
(14) $\begin{split}&d_{{\rm{LoS-
NLoS}},j,k}\left(r\right)=\max\left\\{{{h_{k}},{{\left({\frac{{{\eta_{{\rm{NLoS}}}}{\rho_{k}}\mathbb{E}[{G_{{\rm{NLoS}}}}]}}{{{\eta_{{\rm{LoS}}}}{\rho_{j}}\mathbb{E}[{G_{{\rm{LoS}}}}]}}}\right)}^{\frac{1}{\alpha_{{\rm{NLoS}}}}}}{r^{\frac{\alpha_{{\rm{LoS}}}}{\alpha_{\rm{NLoS}}}}}}\right\\}.\end{split}$
(15)

###### Proof.

According to the received power in (7), the average received power of the
associated UAV at distance$r$in tier $j$ is
$S_{j}={\eta_{\rm{LoS}}}{\rho_{j}}{G_{{\rm{LoS}}}}{r^{-{\alpha_{{\rm{LoS}}}}}}$.
The closest interfering LoS UAV in tier $k$ is at least at distances
$d_{{\rm{LoS-LoS}},j,k}$, which can be obtained by solving the equality
${\eta_{\rm{LoS}}}{\rho_{j}}{G_{{\rm{LoS}}}}{r^{-{\alpha_{{\rm{LoS}}}}}}={\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}d_{{\rm{LoS-
LoS}},j,k}^{-{\alpha_{{\rm{LoS}}}}}$. Similarly, for NLoS interfering UAVs,
$d_{{\rm{LoS-NLoS}},j,k}$ can be obtained by solving the equality
${\eta_{\rm{LoS}}}{\rho_{j}}{G_{{\rm{LoS}}}}{r^{-{\alpha_{{\rm{LoS}}}}}}={\eta_{{\rm{NLoS}}}}{\rho_{k}}{G_{{\rm{NLoS}}}}d_{{\rm{LoS-
NLoS}},j,k}^{-{\alpha_{{\rm{NLoS}}}}}$.

∎

######  Lemma 2.

Given that the typical user is associated with a NLoS UAV located at distance
$r$ in tier $j$, the closest interfering LoS and NLoS UAV in tier $k$ are at
least at distances ${d_{{\rm{NLoS}}{\rm{-}}{\rm{LoS}},j,k}}\left(r\right)$ and
${d_{{\rm{NLoS}}{\rm{-}}\rm{NLoS},j,k}}\left(r\right)$ , given by

$\begin{split}&{d_{{\rm{NLoS-
LoS}},j,k}}\left(r\right)=\max\left\\{{{h_{k}},{{\left({\frac{{{\eta_{{\rm{LoS}}}}{\rho_{k}}\mathbb{E}[{G_{{\rm{LoS}}}}]}}{{{\eta_{{\rm{NLoS}}}}{\rho_{j}}\mathbb{E}[{G_{{\rm{NLoS}}}}]}}}\right)}^{\frac{1}{{{\alpha_{{\rm{LoS}}}}}}}}{r^{\frac{{{\alpha_{{\rm{NLoS}}}}}}{{{\alpha_{{\rm{LoS}}}}}}}}}\right\\},\end{split}$
(16) ${d_{{\rm{NLoS-
NLoS}},j,k}}\left(r\right)=\max\left\\{{{h_{k}},{{\left({\frac{{{\rho_{k}}}}{{{\rho_{j}}}}}\right)}^{\frac{1}{{{\alpha_{{\rm{NLoS}}}}}}}}r}\right\\}.$
(17)

###### Proof.

The proof is similar to that of lemma 1. ∎

In the subsequent analysis, the horizontal distance is more practical than the
Euclidean distance in this model. The horizontal distances
${z_{Q,j,k}}\left(r\right)$ corresponding to lemma 1 and lemma 2 are defined
as

${z_{Q,j,k}}\left(r\right)=\sqrt{d_{{}_{Q,j,k}}^{2}\left(r\right)-h_{k}^{2}},$
(18)

where
$Q=\\{{\rm{LoS}}-{\rm{LoS}},{\rm{LoS}}-{\rm{NLoS}},{\rm{NLoS}}-{\rm{LoS}},{\rm{NLoS}}-{\rm{NLoS}}\\}$.

###  III-B Distance Distribution of Tagged UAV

Before deriving the PDF of the distance of associated UAV, obtaining the
distance distributions of tagged UAVs is necessary. The distance distributions
are given in the following lemmas.

######  Lemma 3.

Given the distance between the typical user and the origin is $z_{u}$, the CDF
of the distance between the tagged LoS UAV in tier $k$ and the typical user is
given by,

$\small\begin{split}&{F_{{\rm{RLoS}},k}}\left({r,z_{u}}\right)=1-\exp\Bigg{(}-\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\int_{-\varphi_{\rm{LoS-
LoS}}\left({l,r,{z_{u}}}\right)}^{\varphi_{\rm{LoS-
LoS}}\left({l,r,{z_{u}}}\right)}{v_{k}^{\rm{LoS}}\left({z_{u}},l,\theta\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)},\end{split}$
(19)

where

$\varphi_{Q,j,k}\left({l,r,{z_{u}}}\right)=\arccos\left({\frac{{{l^{2}}+z_{u}^{2}-z_{Q,j,k}^{2}\left(r\right)}}{{2\,l\,{z_{u}}}}}\right),$
(20)
$v_{k}^{Q}\left({z_{u}},l,\theta\right)=\left|l\right|\,{\Lambda_{\rm{UAV},k}}\left(l\right)\,P_{k}^{Q}\left(d_{\rm{u2U}}\left({z_{u}},l,\theta\right)\right),$
(21)

where
$Q=\\{{\rm{LoS}}-{\rm{LoS}},{\rm{LoS}}-{\rm{NLoS}},{\rm{NLoS}}-{\rm{LoS}},{\rm{NLoS}}-{\rm{NLoS}}\\}$,
the horizontal distances ${z_{Q,j,k}}\left(r\right)$ are defined in (18),
$\Lambda_{\rm{UAV},k}\left(r\right)$ and ${P_{k}^{{\rm{LoS}}}\left(z\right)}$
are given in (2) and (3), respectively. The distance between the potential
interfering UAV and the typical user
${d_{\rm{u2U}}}\left({{z_{u}},l,\theta}\right)$ in (21) is given by,

${d_{\rm{u2U}}}\left({{z_{u}},l,\theta}\right)=\sqrt{{{\left({{z_{u}}-l\cos\theta}\right)}^{2}}+{{\left({l\sin\theta}\right)}^{2}}}.$
(22)

###### Proof.

See Appendix A. ∎

######  Lemma 4.

Given the distance between the typical user and the origin is $z_{u}$, the CDF
of distance between the tagged NLoS UAV in tier $k$ and the typical user is
given by

$\small\begin{split}&{F_{{\rm{RNLoS}},k}}\left({r,z_{u}}\right)=1-\exp\Bigg{(}-\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\int_{-\varphi_{\rm{NLoS-
NLoS}}\left({l,r,{z_{u}}}\right)}^{\varphi_{\rm{NLoS-
NLoS}}\left({l,r,{z_{u}}}\right)}{v_{k}^{\rm{NLoS}}\left({z_{u}},l,\theta\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)},\end{split}$
(23)

where ${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ and
$\varphi_{Q,j,k}\left({l,r,{z_{u}}}\right)$ are given in (21) and (20),
respectively.

###### Proof.

The proof is similar to that of Lemma 3, therefore omitted here. ∎

######  Lemma 5.

Given the distance between the typical user and the origin is $z_{u}$, the PDF
of distance between the tagged $Q$ UAV in tier $k$ and the typical user is
given by,

$\begin{split}&{f_{{\rm{R}Q},k}}\left({r,z_{u}}\right)=\exp\Bigg{(}-\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\int_{-\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}v_{k}^{Q}\left({z_{u}},l,\theta\right){\mathsf{d}\theta\mathsf{d}l}\Bigg{)}\\\
&\times\Bigg{(}\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}-\frac{4r\,\mathbbm{1}\left({r>{h_{k}}}\right)\,v_{k}^{Q}\left({z_{u}},l,\varphi_{Q-Q}\right)}{\sqrt{4\,{l^{2}}\,z_{u}^{2}-{{\left({{l^{2}}+z_{u}^{2}-{r^{2}}+h_{k}^{2}}\right)}^{2}}}}\mathsf{d}l\\\
&+\int_{-\varphi_{Q-Q}\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}},r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}},r,{z_{u}}}\right)}\frac{r\;v_{k}^{Q}\left({z_{u}},{z_{u}+\sqrt{{r^{2}}-h_{k}^{2}}},\theta\right)}{\sqrt{{r^{2}}-h_{k}^{2}}}\mathsf{d}\theta\\\
&+\int_{-\varphi_{Q-Q}\left({{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}},r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}},r,{z_{u}}}\right)}\frac{r\;v_{k}^{Q}\left({z_{u}},{z_{u}-\sqrt{{r^{2}}-h_{k}^{2}}},\theta\right)}{\sqrt{{r^{2}}-h_{k}^{2}}}\mathsf{d}\theta\Bigg{)},\\\
\end{split}$ (24)

where $\mathbbm{1}\left({r>{h_{k}}}\right)$ is an indicator function, its
value is $1$ when $r>{h_{k}}$ is satisfied, otherwise $0$,
${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ and
$\varphi_{Q,j,k}\left({l,r,{z_{u}}}\right)$ are given in (21) and (20),
respectively. For a LoS tagged UAV, $Q$ in
(LABEL:sequation:_PDF_of_tagged_UAV) is replaced with LoS, while $Q$ is
replaced with NLoS for an NLoS tagged UAV.

###### Proof.

See Appendix B. ∎

###  III-C Association Probabilities

Association probability is used to describe the probability that a tagged UAV
will eventually be selected as the associated UAV. For the tagged LoS UAV in
tier $k$, there will be no LoS UAVs providing stronger power in tier $k$ than
the tagged UAV, while the NLoS UAVs in tier $k$ may provide stronger average
received power, and in other tiers, both LoS and NLoS UAVs may provide
stronger power. As a result, the association probabilities are given in the
following lemmas.

######  Lemma 6.

Given the distance between the typical user and the origin is $z_{u}$, for the
LoS tagged UAV from tier $j$ at Euclidean distance $r$ from the typical user,
the probability that the typical user is associated with this specific UAV is
given by

$\begin{split}P_{{\rm{LoS}},j}^{A}\left({r,z_{u}}\right)&=\prod\limits_{k=1,j\neq
k}^{K}\exp\Bigg{(}-\int_{{z_{u}}-{z_{{\rm{LoS-
LoS}},j,k}}\left(r\right)}^{{z_{u}}+{z_{{\rm{LoS-
LoS}},j,k}}\left(r\right)}\int_{-\varphi_{\rm{LoS-
LoS}}\left({l,r,{z_{u}}}\right)}^{\varphi_{\rm{LoS-
LoS}}\left({l,r,{z_{u}}}\right)}{v_{k}^{{\rm{LoS}}}\left({{z_{u}},l,\theta}\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)}\\\
&\times\prod\limits_{k=1}^{K}\exp\Bigg{(}-\int_{{z_{u}}-{z_{{\rm{LoS-
NLoS}},j,k}}\left(r\right)}^{{z_{u}}+{z_{{\rm{LoS-
NLoS}},j,k}}\left(r\right)}\int_{-\varphi_{\rm{LoS-
NLoS}}\left({l,r,{z_{u}}}\right)}^{\varphi_{\rm{LoS-
NLoS}}\left({l,r,{z_{u}}}\right)}{v_{k}^{{\rm{NLoS}}}\left({{z_{u}},l,\theta}\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)},\end{split}$
(25)

where ${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ and
$\varphi_{Q,j,k}\left({l,r,{z_{u}}}\right)$ are given in (21) and (20),
respectively.

###### Proof.

See Appendix C. ∎

######  Lemma 7.

Given the distance between the typical user and the origin is $z_{u}$, for the
NLoS tagged UAV from tier $j$ at Euclidean distance $r$ from the typical user,
the probability that the typical user is associated with this specific UAV is
given by

$\begin{split}P_{{\rm{NLoS}},j}^{A}&\left({r,z_{u}}\right)=\prod\limits_{k=1}^{K}\exp\Bigg{(}-\int_{{z_{u}}-{z_{{\rm{NLoS-
LoS}},j,k}}\left(r\right)}^{{z_{u}}+{z_{{\rm{NLoS-
LoS}},j,k}}\left(r\right)}\int_{-\varphi_{\rm{NLoS-
LoS}}\left({l,r,{z_{u}}}\right)}^{\varphi_{\rm{NLoS-
LoS}}\left({l,r,{z_{u}}}\right)}{v_{k}^{{\rm{LoS}}}\left({z_{u}},l,\theta\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)}\\\
&\times\prod\limits_{k=1,j\neq k}^{K}\exp\Bigg{(}-\int_{{z_{u}}-{z_{{\rm{NLoS-
NLoS}},j,k}}\left(r\right)}^{{z_{u}}+{z_{{\rm{NLoS-
NLoS}},j,k}}\left(r\right)}\int_{-\varphi_{\rm{NLoS-
NLoS}}\left({l,r,{z_{u}}}\right)}^{\varphi_{\rm{NLoS-
NLoS}}\left({l,r,{z_{u}}}\right)}{v_{k}^{{\rm{NLoS}}}\left({z_{u}},l,\theta\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)},\end{split}$
(26)

where ${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ and
$\varphi_{Q,j,k}\left({l,r,{z_{u}}}\right)$ are given in (21) and (20),
respectively.

###### Proof.

The proof is similar to that of Lemma 6, therefore omitted here. ∎

###  III-D Coverage Probability

As an indispensable intermediate result to enable computing coverage
probability, the Laplace Transform of interference is given in the following
lemma.

######  Lemma 8.

Given that the distance between the typical user and the origin is $z_{u}$,
the Laplace transform of the interference power conditioned on the associated
UAV in tier $j$ with Euclidean distance $r$ from the typical user is given by,

$\begin{split}&{\mathcal{L}_{I_{Q_{1},j}}}\left({s,r,z_{u}}\right)=\prod\limits_{k=1}^{K}\Big{[}{\mathcal{L}_{I_{Q_{1}-\rm{LoS},j,k}}}\left({s,r,z_{u}}\right)\times{\mathcal{L}_{I_{Q_{1}-\rm{NLoS},j,k}}}\left({s,r,z_{u}}\right)\Big{]},\end{split}$
(27)

where $Q_{1}$ is replaced with LoS when the typical user is associated with
LoS UAV, $Q_{1}$ is replaced with NLoS when NLoS UAV is associated, and
${\mathcal{L}_{I_{Q_{1}-Q_{2},j,k}}}\left({s,r,z_{u}}\right)$,
$Q_{2}=\\{\rm{LoS},\rm{NLoS}\\}$ is given by,

$\begin{split}{\mathcal{L}_{I_{Q_{1}-Q_{2},j,k}}}\left({s,r,z_{u}}\right)=\exp&\Bigg{(}-\int_{0}^{\max\left\\{{0,{z_{u}}-{z_{Q_{1}-Q_{2},j,k}}\left(r\right)}\right\\}}\int_{-\pi}^{\pi}{v_{k}^{Q_{2}}\left({{z_{u}},l,\theta}\right){w_{Q_{2},k}}\left({s,{z_{u}},l,\theta}\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)}\\\
&\times\exp\Bigg{(}-\int_{{z_{u}}+{z_{Q_{1}-Q_{2},j,k}}\left(r\right)}^{+\infty}\int_{-\pi}^{\pi}{v_{k}^{Q_{2}}\left({{z_{u}},l,\theta}\right){w_{Q_{2},k}}\left({s,{z_{u}},l,\theta}\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)}\\\
\times\exp\Bigg{(}-2&\int_{{z_{u}}-{z_{Q_{1}-Q_{2},j,k}}\left(r\right)}^{{z_{u}}+{z_{Q_{1}-Q_{2},j,k}}\left(r\right)}\int_{\varphi_{Q_{1}-Q_{2},j,k}\left({l,r,{z_{u}}}\right)}^{\pi}{v_{k}^{Q_{2}}\left({{z_{u}},l,\theta}\right){w_{Q_{2},k}}\left({s,{z_{u}},l,\theta}\right)\\!\mathsf{d}\theta\mathsf{d}l}\Bigg{)},\end{split}$
(28)

where

$\begin{split}&{w_{Q_{2},k}}\left({s,r}\right)=1-{\left({\frac{m_{Q_{2}}}{{m_{Q_{2}}}+s{\eta_{Q_{2}}}{\rho_{k}}{G_{Q_{2}}}{{\left({d_{u2U}^{2}\left({{z_{u}},l,\theta}\right)+h_{k}^{2}}\right)}^{\frac{-\alpha_{Q_{2}}}{2}}}}}\right)^{\frac{m_{Q_{2}}}{2}}},\end{split}$
(29)

${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ and
$\varphi_{Q,j,k}\left({l,r,{z_{u}}}\right)$ are given in (21) and (20),
respectively. lemma that $Q_{1}$ represents the type of associated UAV, while
$Q_{2}$ represents the type of interfering UAVs.

###### Proof.

See Appendix D. ∎

As the distance distributions of tagged UAVs and the association probabilities
have been derived, we are ready to calculate the local coverage probability.
The definition and derivation of the local coverage probability are given as
follows.

######  Definition 1 (Local coverage probability).

The local coverage probability ${P^{C}}\left({z_{u},\gamma}\right)$ is the
probability that the SINR of the typical user at distance $z_{u}$ from the
origin is greater than threshold $\gamma$.

######  Theorem 1.

The exact coverage probability ${P^{C}}\left({z_{u},\gamma}\right)$ for the
typical user is given by,

$\small\begin{split}{P^{C}}\left({z_{u},\gamma}\right)&=\sum\limits_{k=1}^{K}\int_{{h_{k}}}^{+\infty}{{f_{{\rm{RLoS}},k}}\left({r,z_{u}}\right)}P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)\sum\limits_{n=0}^{{m_{{\rm{LoS}}}}-1}{\left[{\frac{{{{\left({-s}\right)}^{n}}}}{{n!}}\frac{{{\partial^{n}}}}{{\partial{s^{n}}}}{\mathcal{L}_{U_{\rm{LoS},k}}}\left({s,r,z_{u}}\right)}\right]}_{s={\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)}\mathsf{d}r\\\
&+\sum\limits_{k=1}^{K}\int_{{h_{k}}}^{+\infty}{f_{{\rm{RNLoS}},k}}\left({r,z_{u}}\right)P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)\sum\limits_{n=0}^{{m_{{\rm{NLoS}}}}-1}{{{\left[{\frac{{{{\left({-s}\right)}^{n}}}}{{n!}}\frac{{{\partial^{n}}}}{{\partial{s^{n}}}}{\mathcal{L}_{U_{\rm{NLoS},k}}}\left({s,r,z_{u}}\right)}\right]}_{s={\mu_{{\rm{NLoS}}}}\left({r,\gamma}\right)}}}\mathsf{d}r,\end{split}$
(30)
${\mathcal{L}_{U_{Q,k}}}\left(s,r,z_{u}\right)=\exp\left({-{\sigma^{2}}s}\right){\mathcal{L}_{I_{Q,k}}}\left(s,r,z_{u}\right),$
(31)
${\mu_{Q,k}}\left(r\right){=}m_{Q}\gamma\eta_{Q}^{-1}\rho_{k}^{-1}{r^{{\alpha_{Q}}}},$
(32)

$Q=\\{{\rm{LoS}},{\rm{NLoS}}\\}$ in (31) and (32),
${f_{{\rm{R}}Q,k}}\left({r,z_{u}}\right)$,
${P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)}$ and
${P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)}$ are defined in
(LABEL:sequation:_PDF_of_tagged_UAV), (25) and (26).

###### Proof.

See Appendix E. ∎

Based on the local coverage probability, the definition and derivation of the
overall coverage probability are given as follows.

######  Definition 2 (Overall coverage probability).

The overall coverage probability is the average coverage probability of all
users.

From the definition, the overall coverage probability for the typical user is
the normalized expectation of the local coverage probability with regard to
$z_{u}$.

######  Corollary 1.

The overall exact coverage probability with the SINR threshold $\gamma$ is
given by,

${P_{\rm{Overall}}^{C}}\left(\gamma\right)=\frac{{\int_{0}^{+\infty}{\Lambda_{u}}\left({{z_{u}}}\right){{P^{C}}\left({{z_{u}},\gamma}\right){z_{u}}\mathsf{d}{z_{u}}}}}{{\int_{0}^{+\infty}{{\Lambda_{u}}\left({{z_{u}}}\right){z_{u}}\mathsf{d}{z_{u}}}}}.$
(33)

As is shown in (30), higher-order derivatives of the Laplace transform are
needed while deriving the exact coverage probability. Because the
computational complexity increases rapidly as the order of the derivative
increases, the amount of computation is not acceptable under large shape
parameters $m_{\rm{LoS}}$ and $m_{\rm{NLoS}}$. Therefore, we provide an
approximate evaluation of the coverage probability using the upper bound of
the CDF of the Gamma distribution [28].

######  Theorem 2.

The approximate coverage probability
${\widetilde{P}^{C}}\left({z_{u},\gamma}\right)$ for the typical user is given
by,

$\small\begin{split}{\widetilde{P}^{C}}&\left({z_{u}},\gamma\right)=\sum\limits_{k=1}^{K}\int_{{h_{k}}}^{+\infty}{{f_{{\rm{RLoS}},k}}\left({r,z_{u}}\right)}P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)\sum\limits_{n=1}^{{m_{{\rm{LoS}}}}}{\binom{m_{{\rm{LoS}}}}{n}{{\\!\left({-1}\right)}^{n+1}}\\!\\!{\mathcal{L}_{U_{\rm{LoS},k}}}\\!\\!\left({n\,{\omega_{\rm{LoS}}}\,{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)},r,z_{u}\right)}\mathsf{d}r\\\
&\\!\\!+\sum\limits_{k=1}^{K}\int_{{h_{k}}}^{+\infty}{f_{{\rm{RNLoS}},k}}\left({r,z_{u}}\right)P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)\sum\limits_{n=1}^{{m_{{\rm{NLoS}}}}}\\!\\!{\binom{m_{{\rm{NLoS}}}}{n}{{\\!\left(\\!{-1}\\!\right)}^{n+1}}\\!\\!{\mathcal{L}_{U_{\rm{NLoS},k}}}\\!\\!\left({n\,{\omega_{\rm{NLoS}}}\,{\mu_{{\rm{NLoS}}}}\\!\left({r,\gamma}\right)},\\!r,\\!z_{u}\right)}\mathsf{d}r,\end{split}$
(34)

where

$\omega_{Q}=\left({m_{Q}}!\right)^{-\frac{1}{m_{Q}}},Q=\\{{\rm{LoS}},{\rm{NLoS}}\\},$
(35)

${f_{{\rm{R}}Q,k}}\left({r,z_{u}}\right)$,
${P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)}$,
${P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)}$ and ${\mu_{Q,k}}$ are defined
in (LABEL:sequation:_PDF_of_tagged_UAV), (25), (26) and (32).

###### Proof.

See Appendix F. ∎

The same as the overall exact coverage probability, the overall approximate
coverage probability is given in the following corollary.

######  Corollary 2.

The overall approximate coverage probability is given by,

${\widetilde{P}_{\rm{Overall}}^{C}}\left(\gamma\right)=\frac{{\int_{0}^{+\infty}{\Lambda_{u}}\left({{z_{u}}}\right){{\widetilde{P}^{C}}\left({{z_{u}},\gamma}\right){z_{u}}\mathsf{d}{z_{u}}}}}{{\int_{0}^{+\infty}{{\Lambda_{u}}\left({{z_{u}}}\right){z_{u}}\mathsf{d}{z_{u}}}}}.$
(36)

##  IV Numerical Results

In this section, we compare the coverage performance of different systems and
optimize the overall coverage probability by changing the distribution of UAVs
in different tiers. Referring to [9, 17, 29], we assume the channel parameters
as follows: the LoS and NLoS path-loss exponents are $\alpha_{\rm{LoS}}=2$ and
$\alpha_{\rm{NLoS}}=3$, the mean additional gains for LoS and NLoS
transmissions are $\eta_{\rm{LoS}}=0\rm{dB}$ and
$\eta_{\rm{NLoS}}=-20\rm{dB}$, $m$ parameters of Nakagami-m fading for LoS and
NLoS UAVs are $m_{\rm{LoS}}=2$ and $m_{\rm{NLoS}}=1$, the noise power is
$\sigma^{2}=10^{-7}\rm{W}$, the parameters for the probability of establishing
an LoS link in (3) are $a=4.88$ and $b=0.429$. The deterministic parameters of
users’ distribution in (1) are $\lambda_{u}=10^{-3}\rm{m}^{-2}$ and
$\beta_{u}=5\times 10^{-3}$. As a non-homogeneous PPP, the distribution of
users can be realized by thinning property [11] of homogeneous PPP. Finally,
we assume three tiers of UAVs are deployed in a small town center square with
sides of 5km, at 50,100 and 150 meters height, with the corresponding
transmission power 2, 7, and 12 dBm, respectively, and the same value of
$\lambda_{1}=\lambda_{2}=\lambda_{3}=4\times 10^{-5}\rm{m}^{-2}$.

TABLE I: Table of Parameters |  $h_{k}$ [m] | $\beta$ | $P^{C}_{\tiny{\rm{Overall}}}$ | Density of UAVs
---|---|---|---|---
One-tier | 50 | 3.2 $\times 10^{-3}$ | 0.9026 | $1\,\rm{UAV}/\rm{km}^{2}$
One-tier | 100 | 3.2 $\times 10^{-3}$ | 0.9203 | $1\,\rm{UAV}/\rm{km}^{2}$
One-tier | 150 | 3.2 $\times 10^{-3}$ | 0.9367 | $1\,\rm{UAV}/\rm{km}^{2}$
Three-tier | (50,100,150) |  $(4.5,5.8,7.6)$ $\times 10^{-3}$ | 0.9713 | $1\,\rm{UAV}/\rm{km}^{2}$
Uniform Distribution | 100 | 0 | 0.3845 | $1\,\rm{UAV}/\rm{km}^{2}$

In Fig. 2, a curve of local coverage probability for the typical user as a
function of the distance between the origin and the typical user is plotted.
We compare the coverage performance of population density-inspired UAV systems
(one-tier and three-tier) with that of the uniformly distributed UAV system.
All of the above systems have the same UAV deployment density on average as
$25\,\rm{UAVs}/km^{2}$ (i.e., $\lambda_{h}=10^{-6}\rm{m}^{-2}$ and
$\beta_{h}=0$ for uniform distribution). The distribution parameters $\beta$
of the other three systems are shown in the table, $\lambda=4\times
10^{-5}\rm{m}^{-2}$ as mentioned above.

Figure 2: Local Coverage Probability for the Typical User.

As shown in Table I, different $\beta$ are chosen to keep the density of the
UAVs as $1\,\rm{UAV}/\rm{km}^{2}$. Assume the threshold of coverage
probability is $\gamma=-15\rm{dB}$. It can be seen that the results of the
Monte-Carlo simulation (lines) coincide well with the results of the
theoretical analysis (points). With the same number of UAVs, no matter how far
away the typical user is from the origin, the performance of the three-tier
network is always better than that of the one-tier networks. As shown in Table
I, the three-tier network has significant advantages in terms of overall
coverage. At the edge, the advantage of the three-tier network is further
expanded, indicating that the lower limit of network coverage or SINR can be
guaranteed. In addition, the coverage performance of all three systems is
declining from center to edge due to the reduced density of UAV deployment.
Compared with the proposed distribution, the uniform distribution system only
has a slight advantage in the edge area, but the overall performance is far
inferior to the other four systems. The clustering effect brings the non-
uniform distribution advantages over the uniform distribution. According to
Fig. 2, when the user is close to the center, the coverage probability of the
UAV network under resident population density-inspired distribution is
significantly greater than that under the uniform distribution. Most users are
clustered in the area close to the center, so the proposed distribution has
significant advantages in the overall coverage probability.

We study the following optimization problems and record the results in Fig. 3
and Fig. 4. We want to maximize the overall coverage probability by changing
the distribution of UAVs, under the premise that the total number of UAVs is
limited (the first constraint) and the coverage of users in any location is
guaranteed (the second constraint). Therefore, the mathematical representation
of the optimization problem is as follows,

$\displaystyle\underset{\beta_{1},\beta_{2}}{\arg\max}$ $\displaystyle
P^{C}_{{\rm Overall}}\left(\gamma_{1}\right)$ (37) $\displaystyle{\rm s.t.}$
$\displaystyle\sum_{k=1}^{K}\int_{0}^{+\infty}{\Lambda_{{\rm
UAV},k}}\left({z}\right)\mathrm{d}z\leq N_{\max},$
$\displaystyle\mathbb{P}\left[\rm{SINR}\geq\gamma_{2}\right]\geq 0.95,\ \ \ \
\forall z_{u}\geq 0,$

where the threshold of overall coverage probability is $\gamma_{1}=-8\rm{dB}$
and the threshold of local coverage for the typical user is
$\gamma_{2}=-20\rm{dB}$. The total number of UAVs is limited to
$N_{\max}=1000$, that is, the maximum density of UAV is
$40\,\rm{UAVs}/\rm{km}^{2}$. We deploy UAVs at the first the second tiers
($h=50,100\,\rm{m}$).

Figure 3: Overall coverage probability under different UAV distributions.

A paradoxical but interesting conclusion in Fig. 3 is, the UAV distribution
preferred by the system is a non-uniform one influenced by the resident
density distribution, but the optimal distribution is not closely related to
the residents numerically. The optimal overall coverage probability
$P_{\small\rm{Overall}}^{C*}$ and the corresponding $\beta_{1}^{*}$ and
$\beta_{2}^{*}$ value are marked in the figure. It can be seen that there is a
considerable difference among $\beta_{1}^{*}$, $\beta_{2}^{*}$ and
$\beta_{u}$. For a large $\beta$, the first constraint cannot be satisfied due
to the high density of UAVs. Under this condition, the coverage probability is
set to 0, so the dark blue area at the bottom left appears. When $\beta$
increases to $10^{-3}$, there is a large amount of interference near the
centre because the system still tends to be uniformly distributed and the
density is relatively high. With the increase of $\beta$, the overall coverage
probability is improved rapidly. Near the optimal area, the coverage
performance of the system is no longer very sensitive to both $\beta_{1}$ and
$\beta_{2}$. This is interesting because in such a tolerant system, we do not
have to select an accurate set of optimal parameters to determine the
distribution of UAVs, but only to estimate a range. For a large $\beta$, the
small number of drones are almost all concentrated in the central area, making
it difficult for users in the edge area to maintain good communication
conditions. Therefore, the second constraint cannot be satisfied, and the dark
blue area appears at the top right of the image. Finally, with the same number
of UAVs as the optimal distribution, the overall coverage probability of the
uniform distribution in the same condition ($h=50,100\,\rm{m}$,
$\lambda_{1}=\lambda_{2}=4\times 10^{-6}\,\rm{m}^{-2}$) is only 0.2883, which
is much lower than 0.8503.

Figure 4: The influence of the UAV distribution in single tier on the overall
coverage probability.

Although the optimization problem (37) has been carefully studied in Fig. 3,
it is still necessary to study the behavior of the system hidden in the dark
blue area. Fig. 4 shows the influence of the UAV distribution in a single tier
on the overall coverage probability. We observe the special case of
$\beta_{2}=10^{-2}$ in Fig. 3 and broaden the range of $\beta_{1}$.

First, it is easy to find that the overall coverage probability increases at
the beginning and then decreases with the increase of the value of
$\beta_{1}$. This can simply be explained by the fact that too many UAVs will
cause too much interference in the central area, while having too few UAVs
will make it difficult for the user to find a close UAV to establish an LoS
link. When $\beta_{1}\leq 10^{-5}$, the number of UAVs in the first tier is
much larger than that in the second tier, so the overall coverage probability
is stable and tends to be similar to that of uniform distribution in the first
tier. When $\beta_{1}\geq 5.6\times 10^{-2}$, the number of UAVs in the first
tier is rapidly decreasing, which means there are no available LoS UAVs nearby
for some users. It is not hard to predict that the final result will converge
to the scenario where only the second tier of UAVs are providing the service.

##  V Further Remarks

###  V-A Analytic Framework Extension

This subsection presents how to extend the existing analysis framework to
other scenarios and network models. Enhancing the coverage is one of the
application scenarios for UAV networks. UAV networks can also relieve the
pressure of insufficient channel capacity in town centers.

######  Remark 1.

Based on Shannon’s theorem and the definition of coverage probability in (13),
the channel capacity can be expressed as [30],

$\mathbb{P}\left[B\log_{2}\left(1+\mathrm{SINR}\right)>\mathcal{R}\right]=\mathbb{P}\left[\mathrm{SINR}>2^{\frac{\mathcal{R}}{B}}-1\right]=\mathbb{P}\left[\mathrm{SINR}>\widetilde{\gamma}\right],$
(38)

where $\widetilde{\gamma}=2^{\frac{\mathcal{R}}{B}}-1$ is the rate threshold.
By replacing $\widetilde{\gamma}$ into SINR threshold $\gamma$, the local
probability and overall probability that the channel capacity is greater than
the rate threshold can be obtained by
${\widetilde{P}^{C}}\left({z_{u}},\widetilde{\gamma}\right)$ given by (34) and
${\widetilde{P}_{\rm{Overall}}^{C}}\left(\widetilde{\gamma}\right)$ given by
(36), respectively.

Next, we study green communications in a small hot spot area centered on a
base station. Energy efficiency, the number of bits that can be transmitted
per unit of energy consumed, is used as a performance metric for green
communication. For convenience, we calculate the energy efficiency as the
ratio of the number of bits transmitted per unit time (channel capacity) to
the energy consumed per unit time (transmission power). The introduction of
the UAV network allows the central base station to reduce its coverage area,
thereby reducing transmission power and enhancing energy efficiency. The
following remark illustrates how the coverage probability analytic framework
can be applied to the above scenario.

######  Remark 2.

From the definition of energy efficiency, it can be calculated by the ratio of
channel capacity to transmission power,

$\mathbb{P}\left[\frac{B}{\rho_{k}}\log_{2}\left(1+\rm{SINR}\right)>\mathcal{E}\right]=\mathbb{P}\left[\rm{SINR}>2^{\frac{\mathcal{E}\rho_{k}}{B}}-1\right]=\mathbb{P}\left[\rm{SINR}>\widehat{\gamma}\right],$
(39)

By replacing $\widehat{\gamma}=2^{\frac{\mathcal{E}\rho_{k}}{B}}-1$ into SINR
threshold $\gamma$, the local probability and overall probability that the
energy efficiency is greater than the rate threshold can be obtained by
${\widetilde{P}^{C}}\left({z_{u}},\widehat{\gamma}\right)$ given by (34) and
${\widetilde{P}_{\rm{Overall}}^{C}}\left(\widehat{\gamma}\right)$ given by
(36), respectively.

Furthermore, there is no need for dense deployment of UAVs near the base
station in this case. Fortunately, the network model can be easily extended to
the above scenario by adjusting the density distribution of UAVs in (2). Under
the premise that the density (whether of the user or the UAV) is only related
to the distance to the town center, the analytical framework of this paper is
applicable to any distribution model.

###  V-B Distribution Parameter Design

According to the above theorems, it can be seen that the relationship between
the coverage probability and the spatial distribution of UAVs is not
straightforward, and obtaining the optimal parameters by optimization tools is
challenging. Therefore, the following qualitative criteria for parameters
about UAVs’ vertical and horizontal distributions are given. Notice that all
of the remarks have been verified by simulation.

######  Remark 3.

Remarks on altitudes $h_{1},h_{2},\dots,h_{K}$ are given as follow.

  * •

In most cases, UAVs have an optimal altitude, and it is better to deploy the
UAVs near the optimal altitude.

  * •

While facing a low communication quality, the UAVs are suggested to be
distributed at a low altitude so that UAVs are closer to users and users in
the LoS region can be covered.

  * •

In a good communication environment, by increasing the deployment altitude,
more users can establish LoS links with UAVs, therefore, increasing the
coverage probability.

######  Remark 4.

Remarks on the number of tiers $K$ are given as follows.

  * •

With the increase of tiers, more parameters can be optimized so that the
coverage performance can be improved to some extent. A network with fewer
tiers can be considered a special case of a network with more tiers. However,
the improvement in coverage performance is limited when more than three tiers
are applied.

  * •

Consider a more general case where the receivers (users) can be divided into
$M$ classes according to different gain and demodulation capabilities. Multi-
tier distribution has significant advantages over single-tier distribution,
and the number of deployment tiers $K$ is recommended to be larger than $M$.
Assuming that the optimal UAV deployment altitude for the receiver of the
$m$-th class is $h_{m}^{*}$, the height of UAVs is suggested to satisfy
$\min\left\\{h_{1}^{*},h_{2}^{*},\dots,h_{M}^{*}\right\\}<h_{k}<\max\left\\{h_{1}^{*},h_{2}^{*},\dots,h_{M}^{*}\right\\},\forall
k$.

######  Remark 5.

Remarks on homogeneity $\beta$ are given as follows.

  * •

The value of $\beta$ is related to the strength of interference power relative
to noise. When the interference power is significantly stronger than the
environmental noise, it is suggested to choose a smaller $\beta$ to make the
distribution of UAVs more homogeneous and vice versa. Considering that the
value of $\beta_{k}$ will affect the average number of UAVs in hot areas when
$\beta_{k}$ is changed, $\lambda_{k}$ is adjusted to keep the average number
of UAVs unchanged.

  * •

We use the exhaustive search to solve the optimization problem about $\beta$
given in (37), which results in the calculation complexity increases
exponentially with $K$. An improved alternate maximization method can be a
substitution for the exhaustive search as described in [31]. The complexity of
this method is $\mathcal{O}\left(NK^{2}\right)$, where $N$ is the preset
maximum number of rounds. The set of suboptimal parameters is obtained by
optimizing from $\beta_{1},\beta_{2},\dots$ to $\beta_{K}$ in order. When
optimizing $\beta_{k}$, the $\beta_{k}$ is repeatedly reduced by the
predefined step size for at most $N$ times, and one of the $\beta$ in the set
$\\{\beta_{1},\beta_{2},\dots,\beta_{k-1}\\}$ is increased, so that the
coverage probability is maximized when the UAV density is unchanged. The
optimization of $\beta_{k}$ ends when the coverage probability no longer
increases.

TABLE II: Optimization of $K$ and $\beta$. $1\,\rm{UAV}/\rm{km}^{2}$ | One-tier | Two-tier | Three-tier | Five-tier
---|---|---|---|---
$h_{1}=50$m | N/A | N/A | $\beta_{1}=4.5\times 10^{-3}$ | $\beta_{1}=5.4\times 10^{-3}$
$h_{2}=75$m | N/A | N/A | N/A | $\beta_{2}=6.5\times 10^{-3}$
$h_{3}=100$m | N/A | $\beta_{3}=4.2\times 10^{-3}$ | $\beta_{3}=5.8\times 10^{-3}$ | $\beta_{3}=7.8\times 10^{-3}$
$h_{4}=125$m | N/A | N/A | N/A | $\beta_{4}=8.2\times 10^{-3}$
$h_{5}=150$m | $\beta_{5}=3.2\times 10^{-3}$ | $\beta_{5}=5.4\times 10^{-3}$ | $\beta_{5}=7.6\times 10^{-3}$ | $\beta_{5}=9.8\times 10^{-3}$
$P^{C}_{\tiny{\rm{Overall}}}$ | 0.9367 | 0.9557 | 0.9713 | 0.9786

The example in Table II provides further explanation for the above remarks. In
Table II, we compare the coverage performance under different numbers of tiers
$K$. The total density and the density of UAVs in each tier are fixed as
$1\,\rm{UAV}/\rm{km}^{2}$ and $\lambda_{k}=4\times 10^{-5}$. The set of
homogeneity $\\{\beta_{1},\beta_{2},\dots,\beta_{K}\\}$ is obtained by
alternate maximization method. Overall, increasing the number of tiers allows
more parameters to be optimized, thus achieving better coverage performance.
The UAV deployment in the one-iter network ($\beta_{5}=3.2\times 10^{-3}$) can
be regarded as a special case of that of a two-tier network
($\\{\beta_{3},\beta_{5}\\}=\\{+\infty,3.2\times 10^{-3}\\}$), but it is not
optimal. Finally, the gain in coverage probability from deploying more than
three tiers of UAV networks is limited.

##  VI Conclusion and Future Work

In this paper, We studied the coverage performance of multi-tier UAV networks
in a centralized urban model. We first derived the distance distribution of
tagged UAVs and association probability for the selected typical user. Based
on this, the analytical expression of downlink coverage probability is given
and proved to be consistent with the Monte-Carlo simulation results. As a
result, the coverage probability for the typical user and intermediate
products are all related to the distance $z_{u}$. Both the local and total
coverage performance are significantly improved by increasing the number of
UAV network tiers. The urban population density-inspired model has a huge
advantage over the uniform distribution performs. However, too much
concentration of UAVs in the central area will bring more noise to the town
center and fail to maintain communication for users at the edge. Therefore,
how to design the distribution of each tier of UAVs is crucial.

One future research direction is introducing interference and noise mitigation
technologies into the framework based on the proposed resident population
density-inspired model. In urban areas, the relatively dense deployment of
UAVs may cause strong interference. Strong environmental noise in town centers
is also one factor limiting the performance of wireless communication. Under
the SG framework, orthogonal channel [32] and directional antenna gain [33]
can be introduced into the system model respectively to reduce interference
and noise power. In addition, we model the users as a PPP, which means that
the user’s movement is undirected and random. Considering that there is a
directional flow of people in the town [34], analyzing the coverage
probability of the urban system based on SG will be challenging and
application-oriented.

##  Appendix A Proof of Lemma 3

When the distance between the typical user and the origin is fixed, given that
the distance between the tagged LoS UAV in tier $k$ and user
${R_{\rm{LoS},k}}$ is a random valuable, the Cumulative Distribution Function
(CDF) of ${R_{\rm{LoS},k}}$ is given by

$\begin{split}{F_{{\rm{RLoS}},k}}\left({r,z_{u}}\right)&={\mathbb{P}}\left[{{R_{{\rm{LoS}},k}}<r}\right]=1-{\mathbb{P}}\left[{{R_{{\rm{LoS}},k}}>r}\right]\\\
&=1-\mathbb{P}\left[{\mathcal{N}\left({{\mathcal{A}_{k}}\left(r\right)}\right)=0}\right]\overset{(a)}{=}1-\exp\left({-\int_{{\mathcal{A}_{k}}\left(r\right)}{{\Lambda_{\rm{UAV},k}}\left(l\right)l\mathsf{d}l\mathsf{d}\theta}}\right),\end{split}$
(40)

where ${\Lambda_{\rm{UAV},k}}\left(l\right)$ is defined in (2),
$\mathcal{N}\left({{\mathcal{A}_{k}}\left(r\right)}\right)$ in (40) counts the
number of the UAVs in region ${{\mathcal{A}_{k}}\left(r\right)}$, which is a
circle at the height of ${h_{k}}$ centered directly above the typical user
with radius $\sqrt{{r^{2}}-h_{k}^{2}}$, and $(a)$ is given by the property of
the general PPP [35],

$\begin{split}\mathbb{P}\left[{\mathcal{N}\left({{{\mathcal{A}_{k}}\left(r\right)}}\right)=n}\right]=&\exp\left({-\int_{{{\mathcal{A}_{k}}\left(r\right)}}{{\Lambda_{\rm{UAV},k}}\left(l\right)l\,\mathsf{d}l\mathsf{d}\theta}}\right)\frac{{\exp{{\left({-\int_{{{\mathcal{A}_{k}}\left(r\right)}}{{\Lambda_{\rm{UAV},k}}\left(l\right)l\,\mathsf{d}l\mathsf{d}\theta}}\right)}^{n}}}}{{n!}}.\end{split}$
(41)

Figure 5: Vertical Viewed System Schematic Figure.

where $z_{u}$ is the horizontal distance from the typical user to the origin,
$\lambda_{u}$ determines the total density of the plane, $\beta_{u}$ is a
measure of homogeneity.

To integrate formulation in (40) over the region of ${{\mathcal{A}_{k}}}$, the
area is divided into infinite concentric circular arcs centered at the point
which is directly above the origin at the height of $h_{k}$. When the radius
$l$ of the concentric circular arc is fixed, the density function
${\Lambda_{\rm{UAV},k}}$ is a constant, the coordinates of the points on the
arc can be uniquely represented by $\theta$. The bold part of the bottom half
of the Fig. 5 is one of the concentric arcs. The upper bound of $\theta$ can
be obtained from the geometric relations in the Fig. 5, denoted as
$\varphi_{\rm{LoS-LoS}}$, defined in (20), and the lower bound of $\theta$ is
$-\varphi_{\rm{LoS-LoS}}$ because of the symmetry of the circle. Furthermore,
the horizontal distance ${d_{\rm{u2U}}}\left({{z_{u}},l,\theta}\right)$
between the typical user and the point on the arc is defined in (22), which
can also be obtained from simple geometrical relationships. Hence, we have the
following equation

$\begin{split}&\int_{{\mathcal{A}_{k}}\left(r\right)}{{\Lambda_{\rm{UAV},k}}\left(l\right)l\mathsf{d}l\mathsf{d}\theta}=\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\int_{-\varphi_{\rm{LoS-
LoS}}\left({l,r,{z_{u}}}\right)}^{\varphi_{\rm{LoS-
LoS}}\left({l,r,{z_{u}}}\right)}{v_{k}^{{\rm{LoS}}}\left({{z_{u}},l,\theta}\right)\mathsf{d}\theta\mathsf{d}l},\end{split}$
(42)

where ${{\Lambda_{\rm{UAV},k}}\left(l\right)}$ and
${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ are defined in 2) and (21),
respectively. It is important to note that ${l}$ may be negative when the
value of ${z_{{\rm{LoS-LoS}},j,k}}$ is greater than the horizontal distance
${z_{u}}$. Therefore, ${|l|}$ is used in the outer integral.

##  Appendix B Proof of Lemma 5

As in Lemma 5, $Q$ is used to represent the type of tagged UAVs, i.e., $Q$ is
replaced with LoS when an LoS UAV is tagged or NLoS otherwise. By taking the
derivative of ${F_{{\rm{R}Q},k}}\left({r,z_{u}}\right)$, the distribution of
the nearest UAVs in tier $k$ with a distance $r$ from the user is obtained,
which is denoted as ${f_{{\rm{R}Q},k}}\left({r,z_{u}}\right)$,

$\begin{split}&{f_{{\rm{R}Q},k}}\left({r,z_{u}}\right)=\frac{\partial}{{\partial
r}}{F_{{\rm{R}Q},k}}\left({r,z_{u}}\right)\\\ &=\frac{\partial}{{\partial
r}}\Bigg{(}1-\exp\Bigg{(}-\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\int_{-\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}{v_{k}^{Q}\left({{z_{u}},l,\theta}\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)}\Bigg{)}\\\
&\overset{(a)}{=}\exp\Bigg{(}-\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\int_{-\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}{v_{k}^{Q}\left({{z_{u}},l,\theta}\right)\mathsf{d}\theta\mathsf{d}l}\Bigg{)}\\\
&\times\Bigg{(}\int_{{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}^{{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\underbrace{\frac{\partial}{{\partial
r}}f_{\rm{in},k}\left({l,{z_{u}},r,\theta}\right)}_{\small\rm{The\ derivative\
of\ the\ integrand}}\mathsf{d}l\\\
&+\underbrace{\frac{{\partial\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\right)}}{{\partial
r}}f_{\rm{in},k}\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}},{z_{u}},r,\theta}\right)}_{\rm{The\
derivative\ of\ the\ integral\ upper\ bound}}\\\
&-\underbrace{\frac{{\partial\left({{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}}}\right)}}{{\partial
r}}f_{\rm{in},k}\left({{z_{u}}-\sqrt{{r^{2}}-h_{k}^{2}},{z_{u}},r,\theta}\right)}_{\rm{The\
derivative\ of\ the\ integral\ upper\ bound}}\Bigg{)},\end{split}$ (43)

where ${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ and
$\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)$ are defined in (21) and (20),
respectively, and $(a)$ follows Leibnitz’s rule,
$f_{\rm{in},k}\left({l,{z_{u}},r,\theta}\right)$ is the integrand of the outer
integral, given by

$\begin{split}&f_{\rm{in},k}\left({l,{z_{u}},r,\theta}\right)=\int_{-\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}{v_{k}^{Q}\left({{z_{u}},l,\theta}\right)\mathsf{d}\theta}.\end{split}$
(44)

For the derivative of the integral upper bound in (LABEL:fRLoS1),

$\begin{split}&\frac{{\partial\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}}}\right)}}{{\partial
r}}f_{\rm{in},k}\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}},{z_{u}},r,\theta}\right)\\\
&=\frac{{r}}{{\sqrt{{r^{2}}-h_{k}^{2}}}}\int_{-\varphi_{Q-Q}\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}},r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}},r,{z_{u}}}\right)}{v_{k}^{Q}\left({{z_{u}},{z_{u}}+\sqrt{{r^{2}}-h_{k}^{2}},\theta}\right)\mathsf{d}\theta}.\end{split}$
(45)

The derivative of the integral lower bound is similar to that of
(LABEL:fRLoS2), therefore omitted here.

For the derivative of the intergrad,

$\begin{split}&\frac{\partial}{{\partial
r}}f_{\rm{in},k}\left({l,{z_{u}},r,\theta}\right)=\frac{\partial}{{\partial
r}}\int_{-\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}^{\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}{v_{k}^{Q}\left({{z_{u}},l,\theta}\right)\mathsf{d}\theta}\\\
&\overset{(a)}{=}v_{k}^{Q}\left({{z_{u}},l,\varphi_{Q-Q}}\right)\frac{{\partial\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}}{{\partial
r}}-v_{k}^{Q}\left({{z_{u}},l,-\varphi_{Q-Q}}\right)\frac{{\partial\left({-\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}\right)}}{{\partial
r}}\Big{)}\\\
&\overset{(b)}{=}2v_{k}^{Q}\left({{z_{u}},l,\varphi_{Q-Q}}\right)\frac{{\partial\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}}{{\partial
r}}\\\ &\overset{(c)}{=}\mathbbm{1}\left({r>{h_{k}}}\right)\
\frac{{4\,r\,v_{k}^{Q}\left({{z_{u}},l,\varphi_{Q-Q}}\right)}}{{\sqrt{4\,{l^{2}}\,z_{u}^{2}-{{\left({{l^{2}}+z_{u}^{2}-{r^{2}}+h_{k}^{2}}\right)}^{2}}}}},\end{split}$
(46)

where $(a)$ follows Leibniz’s rule for internal integral, and the expression
in $(b)$ is simplified by the fact
${d_{{\rm{u2U}}}}\left({{z_{u}},l,-\varphi_{Q-Q}}\right)={d_{{\rm{u2U}}}}\left({{z_{u}},l,\varphi_{Q-Q}}\right)$,
which can be easily obtained by (20),
$\mathbbm{1}\left({r>{h_{k}}}\right)$ is the indicator function defined in
(5), and $(c)$ is obtained by substitute
$\frac{{\partial\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}}{{\partial r}}$,
which is given by,

$\begin{split}&\frac{{\partial\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)}}{{\partial
r}}=\frac{\partial}{{\partial
u}}\arccos\left({\frac{{{l^{2}}+z_{u}^{2}-{u^{2}}}}{{2\,l\,{z_{u}}}}}\right)\
\frac{\partial}{{\partial v}}\sqrt{{v^{2}}-h_{k}^{2}}\ \frac{{\partial
v}}{{\partial r}}\\\
&=\frac{{2u}}{{\sqrt{4\,{l^{2}}\,z_{u}^{2}-{{\left({{l^{2}}+z_{u}^{2}-{u^{2}}}\right)}^{2}}}}}\frac{{2v}}{{2\sqrt{{v^{2}}-h_{k}^{2}}}}\cdot\mathbbm{1}\left({{{\left({\frac{{{\rho_{k}}}}{{{\rho_{k}}}}}\right)}^{\frac{1}{{{\alpha_{{\rm{LoS}}}}}}}}r>{h_{k}}}\right)\\\
&=\mathbbm{1}\left({r>{h_{k}}}\right)\
\frac{{2r}}{{\sqrt{4\,{l^{2}}\,z_{u}^{2}-{{\left({{l^{2}}+z_{u}^{2}-{r^{2}}+h_{k}^{2}}\right)}^{2}}}}},\end{split}$
(47)

where $u={z_{{\rm{LoS-
LoS}},j,k}}\left(r\right)\left|{{}_{j=k}}\right.=\mathbbm{1}\left({r>{h_{k}}}\right)\,\sqrt{{r^{2}}-h_{k}^{2}}$,
and $v={d_{{\rm{LoS-
LoS}},j,k}}\left(r\right)\left|{{}_{j=k}}\right.=\max\left\\{{{h_{k}},r}\right\\}$.
Substitute (LABEL:fRLoS2) and (46) into (LABEL:fRLoS1), the final result is
derived.

Figure 6: Two relationships between user’s non-interference circle and
interfering UAVs’ circle.

##  Appendix C Proof of Lemma 6

When the distance between the typical user and the origin is fixed, the
probability that the typical user is associated with the tagged LoS UAV in
tier $j$ is equal to the probability that the average received power of other
$2K-1$ tagged UAVs is lower than it, where $K$ is the number of tiers. Using
the solution of lemma 1, we have

$\begin{split}&P_{{\rm{LoS}},j}^{A}({r,z_{u}})={\prod\limits_{k=1,j\neq
k}^{K}{\mathbb{P}}\left[{{R_{{\rm{LoS}},k}}>{d_{{\rm{LoS-
LoS}},j,k}}\left(r\right)}\right]}\times{\prod\limits_{k=1}^{K}{\mathbb{P}}\left[{{R_{{\rm{NLoS},k}}}>{d_{{\rm{LoS-
NLoS}},j,k}}\left(r\right)}\right]}\\\ &={\prod\limits_{k=1,j\neq
k}^{K}{\mathbb{P}}\left[{\mathcal{N}\left({{\mathcal{A}_{\rm{LoS},k}}\left(r\right)}\right)=0}\right]}\times{\prod\limits_{k=1}^{K}{\mathbb{P}}\left[{\mathcal{N}\left({{\mathcal{A}_{\rm{NLoS},k}}\left(r\right)}\right)=0}\right]}\\\
&\overset{(a)}{=}{\prod\limits_{k=1,j\neq
k}^{K}\exp\left({-\int_{{\mathcal{A}_{\rm{LoS},k}}\left(r\right)}{{v_{k}^{\rm{LoS}}\left({{z_{u}},l,\theta}\right)}\mathsf{d}l\,\mathsf{d}\theta}}\right)}\times{\prod\limits_{k=1}^{K}\exp\left({-\int_{{\mathcal{A}_{\rm{NLoS},k}}\left(r\right)}{{v_{k}^{\rm{NLoS}}\left({{z_{u}},l,\theta}\right)}\mathsf{d}l\,\mathsf{d}\theta}}\right)},\end{split}$
(48)

where ${R_{Q,k}}$ is the distance between the tagged UAV in tier $k$ and
typical user, ${v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}$ is defined in (21),
$\mathcal{N}\left({{\mathcal{A}_{Q,k}}\left(r\right)}\right)$ counts the
number of the UAVs in region ${{\mathcal{A}_{Q,k}}\left(r\right)}$, which is a
circle at the height of ${h_{k}}$ centered directly above the typical user
with radius ${z_{\rm{LoS}-Q,j,k}}\left(r\right)$,
$Q=\left\\{\rm{LoS},\rm{NLoS}\right\\}$, and $(b)$ is given by the property of
the general PPP in (41). The following two equations can be obtained in a
similar way to (42),

$\begin{split}&\int_{{\mathcal{A}_{Q,k}}\left(r\right)}{{v_{k}^{Q}\left({{z_{u}},l,\theta}\right)}\mathsf{d}l\,\mathsf{d}\theta}=\int_{{z_{u}}-{z_{{Q-Q},j,k}}}^{{z_{u}}+{z_{Q-Q,j,k}}}\int_{-\varphi{{}_{Q-Q}}\left(l,r,z_{u}\right)}^{\varphi{{}_{Q-Q}}\left(l,r,z_{u}\right)}{v_{k}^{Q}({z_{u},l,\theta}})\mathsf{d}\theta\mathsf{d}l,\end{split}$
(49)

where $Q=\left\\{\rm{LoS},\rm{NLoS}\right\\}$. The final result is derived by
substituting (49) into (48).

##  Appendix D Proof of Lemma 8

For LoS associated UAV in tier $j$, the Laplace transform of the interference
power can be expressed as,

$\begin{split}&{\mathcal{L}_{I_{\rm{LoS},j}}}\left({s,r,z_{u}}\right)\overset{(a)}{=}\mathbb{E}_{I_{\rm{LoS},j}}\left[{{e^{-s{I_{\rm{LoS},j}}}}}\right]\\\
&={\mathbb{E}_{{\Phi},G}}\Bigg{[}\exp\Bigg{(}-s\sum\limits_{k=1}^{K}\Bigg{(}\sum\limits_{x\in{\Phi_{{\rm{LoS}},k}}\backslash{x_{o}}}{{\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}{r^{-{\alpha_{{\rm{LoS}}}}}}}+\sum\limits_{x\in{\Phi_{{\rm{NLoS}},k}}}{{\eta_{{\rm{NLoS}}}}{\rho_{k}}{G_{{\rm{NLoS}}}}{r^{-{\alpha_{{\rm{NLoS}}}}}}}\Bigg{)}\Bigg{)}\Bigg{]}\\\
&\overset{(b)}{=}\prod\limits_{k=1}^{K}{\mathbb{E}_{{\Phi_{{\rm{LoS}},k}}}}\Bigg{[}\prod\limits_{x\in{\Phi_{{\rm{LoS}},k}}\backslash{x_{o}}}{{\mathbb{E}_{{G_{{\rm{LoS}}}}}}\left[{\exp\left({-s\,{\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}{r^{-{\alpha_{{\rm{LoS}}}}}}}\right)}\right]}\Bigg{]}\\\
&\times\prod\limits_{k=1}^{K}{\mathbb{E}_{{\Phi_{{\rm{NLoS}},k}}}}\Bigg{[}\prod\limits_{x\in{\Phi_{{\rm{NLoS}},k}}}{{\mathbb{E}_{{G_{{\rm{NLoS}}}}}}\left[{\exp\left({-s\,{\eta_{{\rm{NLoS}}}}{\rho_{k}}{G_{{\rm{NLoS}}}}{r^{-{\alpha_{{\rm{NLoS}}}}}}}\right)}\right]}\Bigg{]},\end{split}$
(50)

where $(a)$ follows the definition of Laplace transform, $\\\
{{\Phi_{{\rm{LoS}},k}}\backslash{\Phi_{o}}}$ are all of the LoS UAVs in tier
$k$ except for the associated one, and $(b)$ follows the independence of the
point process and the small scale fading, $(c)$ is obtained. For Laplace
transform of the interference caused by LoS UAVs in tier $k$,

$\small\begin{split}\qquad&{{\mathbb{E}_{{\Phi_{{\rm{LoS}},k}}}}\left[{\prod\limits_{x\in{\Phi_{{\rm{LoS}},k}}\backslash{x_{o}}}{{\mathbb{E}_{{G_{{\rm{LoS}}}}}}\left[{\exp\left({-s\,{\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}{r^{-{\alpha_{{\rm{LoS}}}}}}}\right)}\right]}}\right]}\\\
&\overset{(a)}{=}\exp(-\int_{\mathcal{A}_{\rm{LoS},k}^{C}(r)}v_{k}^{\rm{LoS}}\left({z_{u}},l,\theta\right){\left({1-{\mathbb{E}_{{G_{{\rm{LoS}}}}}}\left[{\exp\left({-s\,{\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}{\left({d_{u2U}^{2}\left({{z_{u}},l,\theta}\right)+h_{k}^{2}}\right)^{-{\alpha_{{\rm{LoS}}}}/2}}}\right)}\right]}\right)}\,\mathsf{d}\theta\mathsf{d}l)\\\
&=\exp\Bigg{(}-\int_{\mathcal{A}_{\rm{LoS},k}^{C}(r)}v_{k}^{\rm{LoS}}\left({{z_{u}},l,\theta}\right){\left({1-{{\left({\frac{{{m_{{\rm{LoS}}}}}}{{{m_{{\rm{LoS}}}}+s\,{\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}{\left({d_{u2U}^{2}\left({{z_{u}},l,\theta}\right)+h_{k}^{2}}\right)^{-{\alpha_{\rm{LoS}}}/2}}}}}\right)}^{m_{{\rm{LoS}}}}}}\right)}l\,\mathsf{d}\theta\mathsf{d}l\Bigg{)}\\\
&\overset{(b)}{=}\exp\Bigg{(}-\underbrace{\int_{0}^{\max\left\\{{0,{z_{u}}-{z_{\rm{LoS-
LoS},j,k}}\left(r\right)}\right\\}}\int_{-\pi}^{\pi}{v_{k}^{{\rm{LoS}}}\left({{z_{u}},l,\theta}\right){w_{{{\rm{LoS}}},k}}\left(s,{z_{u}},l,\theta\right)\mathsf{d}\theta\mathsf{d}l}}_{\rm{User^{\prime}s\
non-interference\ circle\ separates\ from\ interfering\ UAVs^{\prime}\
circle}}\Bigg{)}\\\ &\times\exp\Bigg{(}-\underbrace{\int_{{z_{u}}+{z_{\rm{LoS-
LoS},j,k}}\left(r\right)}^{+\infty}\int_{-\pi}^{\pi}{v_{k}^{{\rm{LoS}}}\left({{z_{u}},l,\theta}\right){w_{{\rm{LoS}},k}}\left(s,{z_{u}},l,\theta\right)\mathsf{d}\theta\mathsf{d}l}}_{\rm{User^{\prime}s\
non-interference\ circle\ contained\ by\ interfering\ UAVs^{\prime}\
circle}}\Bigg{)}\\\
&\times\exp\Bigg{(}-\underbrace{2\int_{{{z_{u}}-{z_{\rm{LoS-
LoS},j,k}}\left(r\right)}}^{{z_{u}}+{z_{\rm{LoS-
LoS},j,k}}\left(r\right)}\int_{\varphi_{\rm{LoS-
LoS},j,k}\left({l,r,{z_{u}}}\right)}^{\pi}{v_{k}^{{\rm{LoS}}}\left({{z_{u}},l,\theta}\right){w_{{\rm{LoS}},k}}\left(s,{z_{u}},l,\theta\right)\mathsf{d}\theta\mathsf{d}l}}_{\rm{User^{\prime}s\
non-interference\ circle\ and\ interfering\ UAVs^{\prime}\ circle\
intersect}}\Bigg{)}\\\ \end{split}$ (51)

In equation (51) shown at the top of the next page, $(a)$ follows the PGFL of
inhomogeneous PPP [35], $\mathcal{A}_{\rm{LoS},k}^{C}(r)$ is the complement of
$\mathcal{A}_{\rm{LoS},k}(r)$ in the two dimensional plane at the height of
${h_{k}}$, the definition of $\mathcal{A}_{\rm{LoS},k}(r)$ is described in
Appendix A, ${v_{k}^{Q}\left({z_{u}},l,\theta\right)}$ and
$\varphi_{Q-Q}\left({l,r,{z_{u}}}\right)$ are defined in (21) and (20)
respectively, ${w_{{\rm{LoS,}}k}}\left({s,r}\right)$ defined in (29) is used
to simplify the expression in $(b)$.

As is shown in Fig. 6, there should be non-interfering UAVs inside the green
circle centered at the typical user, the green circle is called the user’s
non-interference circle. The difference between the left and right images is
whether the origin is included by user’s non-interference circle. For a fixed
radius $l$, the circles centered at the origin is used to cover the possible
locations of interfering UAVs with horizontal distance $l$ to the origin,
called the interfering UAV’s circle. These two circles may be separated or
intersected, and sometimes one circle may contain another, shown in step $(b)$
of (51). The Laplace transform of the interference in other conditions is
similar to the process in (51), therefore omitted here.

##  Appendix E Proof of Theorem 1

By the definition of coverage probability in (13), SINR becomes a
deterministic expression only when: (i) the tier where the associated UAV is
located; (ii) LoS or NLoS link constructed by the typical user and associated
UAV; (iii) the distance between the typical user and the origin; (iv) the
Euclidean distance between the typical user and the associated UAV. Therefore,
the coverage probability of the typical user is given by
(LABEL:exact_coverage_probability_1) at the top of next page.

$\begin{split}&{P^{C}}\left({{z_{u}},\gamma}\right)=\sum\limits_{k=1}^{K}{{\mathbb{E}_{r,I}}\left[{P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)\mathbb{P}\left[{\frac{{{\eta_{{\rm{LoS}}}}{\rho_{k}}{G_{{\rm{LoS}}}}{r^{-{\alpha_{{\rm{LoS}}}}}}}}{{I_{\rm{LoS},k}\left({r,z_{u}}\right)+{\sigma^{2}}}}>\gamma}\right]}\right]}\\\
&+\sum\limits_{k=1}^{K}{{\mathbb{E}_{r,I}}\left[{P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)\mathbb{P}\left[{\frac{{{\eta_{{\rm{NLoS}}}}{\rho_{k}}{G_{{\rm{NLoS}}}}{r^{-{\alpha_{{\rm{NLoS}}}}}}}}{{I_{\rm{NLoS},k}\left({r,z_{u}}\right)+{\sigma^{2}}}}>\gamma}\right]}\right]}\\\
&\overset{(a)}{=}\sum\limits_{k=1}^{K}{{\mathbb{E}_{r,U}}\left[{P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)\mathbb{P}\left[{{G_{{\rm{LoS}}}}>{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)U_{\rm{LoS},k}\left({r,z_{u}}\right)}\right]}\right]}\\\
&+\sum\limits_{k=1}^{K}{{\mathbb{E}_{r,U}}\left[{P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)\mathbb{P}\left[{{G_{{\rm{NLoS}}}}>{\mu_{{\rm{NLoS}}}}\left({r,\gamma}\right)U_{\rm{NLoS},k}\left({r,z_{u}}\right)}\right]}\right]}\\\
&\overset{(b)}{=}\sum\limits_{k=1}^{K}\int_{{h_{k}}}^{+\infty}{\mathbb{E}_{U}}\left[{\mathbb{P}\left[{{G_{{\rm{LoS}}}}>{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)U_{\rm{LoS},k}\left({r,z_{u}}\right)}\right]}\right]P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)f_{{\rm{RLoS}},k}\left({r,z_{u}}\right)\mathsf{d}r\\\
&+\sum\limits_{k=1}^{K}\int_{{h_{k}}}^{+\infty}{\mathbb{E}_{U}}\left[{\mathbb{P}\left[{{G_{{\rm{NLoS}}}}>{\mu_{{\rm{NLoS}}}}\left({r,\gamma}\right)U_{\rm{NLoS},k}\left({r,z_{u}}\right)}\right]}\right]P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)f_{{\rm{RNLoS}},k}\left({r,z_{u}}\right)\mathsf{d}r,\end{split}$
(52)

where $P_{{\rm{LoS}},k}^{A}\left({r,z_{u}}\right)$,
$P_{{\rm{NLoS}},k}^{A}\left({r,z_{u}}\right)$ and
$f_{{\rm{RQ}},k}\left({r,z_{u}}\right)$ are given in (25), (26) and
(LABEL:sequation:_PDF_of_tagged_UAV), respectively, $(a)$ is obtained by
substituting
$U_{Q,k}\left({r,z_{u}}\right)=I_{Q,k}\left({r,z_{u}}\right)+{\sigma^{2}}$ and
${\mu_{Q,k}}\left(r,\gamma\right)$ are define in (32) into the former result,
$(b)$ is obtained from the expectation of $r$. In order to get the final
analytical result, the next steps are taken,

$\begin{split}&{\mathbb{E}_{U}}\left[{\mathbb{P}\left[{{G_{{\rm{LoS}}}}>{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right]}\right]\\\
&\overset{(a)}{=}{\mathbb{E}_{U}}\left[{\frac{{{\Gamma_{u}}\left({{m_{{\rm{LoS}}}},{m_{{\rm{LoS}}}}{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right)}}{{\Gamma\left({{m_{{\rm{LoS}}}}}\right)}}}\right]\\\
&\overset{(b)}{=}{\mathbb{E}_{U}}\Big{[}\exp\left({-{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)U\left({r,z_{u}}\right)}\right)\sum\limits_{n=0}^{{m_{{\rm{LoS}}}}-1}{\frac{{{{\left({{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right)}^{n}}}}{{n!}}}\Bigg{]}\\\
&=\sum\limits_{n=0}^{{m_{{\rm{LoS}}}}-1}{\frac{{{{\left({{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)}\right)}^{n}}}}{{n!}}}{\mathbb{E}_{U}}\big{[}\exp\left({-{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right){{\left({{U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right)}^{n}}\big{]}\\\
&\overset{(c)}{=}\sum\limits_{n=0}^{{m_{{\rm{LoS}}}}-1}{{{\left[{\frac{{{{\left({-s}\right)}^{n}}}}{{n!}}\frac{{{\partial^{n}}}}{{\partial{s^{n}}}}{\mathcal{L}_{U_{\rm{LoS},k}}}\left({s,r,z_{u}}\right)}\right]}_{s={\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)}}},\end{split}$
(53)

where $(a)$ follows the complementary cumulative distribution function (CCDF)
of the Gamma distribution
${\overline{F}_{G}}\left(g\right)=\frac{{{\Gamma_{u}}\left({m,mg}\right)}}{{\Gamma\left(m\right)}}$,
where
${\Gamma_{u}}\left({m,mg}\right)=\int_{mg}^{+\infty}{{t^{m-1}}{e^{-t}}dt}$ is
the upper incomplete Gamma function, and $(b)$ follows the definition
$\frac{{{\Gamma_{u}}\left({m,mg}\right)}}{{\Gamma\left(m\right)}}=\exp\left({-g}\right)\sum\limits_{n=0}^{m-1}{\frac{{{g^{n}}}}{{n!}}}$
[36], by the linearity of the expectation operator and

$\begin{split}&{\mathbb{E}_{U}}\left[{\exp\left({-s{U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right){U_{\rm{LoS},k}}{{\left({r,z_{u}}\right)}^{n}}}\right]={\left({-1}\right)^{n}}\frac{{{\partial^{n}}}}{{\partial{s^{n}}}}{\mathcal{L}_{U_{\rm{LoS},k}}}\left({s,r,z_{u}}\right),\end{split}$
(54)

$(c)$ is obtained. The steps of NLoS UAVs are similar to that of LoS UAVs,
therefore omitted here.

##  Appendix F Proof of Theorem 2

Because the first several steps of the proof of approximate coverage
probability are similar to that of exact coverage probability, we start from
formulation (LABEL:exact_coverage_probability_2) step $(a)$,

$\begin{split}&{\mathbb{E}_{U}}\left[{\frac{{{\Gamma_{u}}\left({{m_{{\rm{LoS}}}},{m_{{\rm{LoS}}}}{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r\left|{{r_{u}}}\right.}\right)}\right)}}{{\Gamma\left({{m_{{\rm{LoS}}}}}\right)}}}\right]\\\
&\overset{(a)}{=}1-{\mathbb{E}_{U}}\left[{\frac{{{\Gamma_{l}}\left({{m_{{\rm{LoS}}}},{m_{{\rm{LoS}}}}{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r\left|{{r_{u}}}\right.}\right)}\right)}}{{\Gamma\left({{m_{{\rm{LoS}}}}}\right)}}}\right]\\\
&\overset{(b)}{\approx}1-{\mathbb{E}_{U}}\left[{{{\left({1-\exp\left({-{\beta_{{\rm{LoS}}}}{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right)}\right)}^{{m_{{\rm{LoS}}}}}}}\right]\\\
&\overset{(c)}{=}{\mathbb{E}_{U}}\Bigg{[}\sum\limits_{n=1}^{{m_{{\rm{LoS}}}}}{\binom{m_{\rm{LoS}}}{n}}{{\left({-1}\right)}^{n+1}}\exp\left({-n{\omega_{{\rm{LoS}}}}{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right){U_{\rm{LoS},k}}\left({r,z_{u}}\right)}\right)\Bigg{]}\\\
&=\sum\limits_{n=1}^{{m_{{\rm{LoS}}}}}{\binom{m_{\rm{LoS}}}{n}{{\left({-1}\right)}^{n+1}}{\mathcal{L}_{U_{\rm{LoS},k}}}\left({n{\omega_{{\rm{LoS}}}}{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)}\right)},\end{split}$
(55)

where ${\mathcal{L}_{U_{\rm{LoS},k}}}\left(s,r,k\right)$ is given in (31), and
$s=n{\omega_{{\rm{LoS}}}}{\mu_{{\rm{LoS},k}}}\left({r,\gamma}\right)$,
$\Gamma_{l}\left(m,mg\right)=\int_{0}^{mg}{{t^{m-1}}{e^{-t}}dt}$ in step $(a)$
is the lower incomplete Gamma function, which satisfies
$\frac{{{\Gamma_{u}}\left({m,mg}\right)}}{{\Gamma\left(m\right)}}=1-\frac{{{\Gamma_{l}}\left({m,mg}\right)}}{{\Gamma\left(m\right)}}$.
$(b)$ follows from the tight approximation to coverage probability, where
${\omega_{{\rm{LoS}}}}={\left({{m_{{\rm{LoS}}}}!}\right)^{\frac{{-1}}{{{m_{{\rm{LoS}}}}}}}}$.
It has been proved in [37] that the tighter upper bound provides an accurate
approximation of the CDF of the Gamma distribution, which is bounded by

${\left({1-{e^{-{\omega_{1}}mg}}}\right)^{m}}<\frac{{{\Gamma_{l}}\left({m,mg}\right)}}{{\Gamma\left(m\right)}}<{\left({1-{e^{-{\omega_{2}}mg}}}\right)^{m}},$
(56)

where $m\neq 1$, and

${\omega_{1}}=\left\\{{\begin{array}[]{*{20}{c}}{1,}&{{\rm{if\,\,}}m>1}\\\
{{{\left({m!}\right)}^{\frac{{-1}}{m}}},}&{{\rm{if\,\,}}m<1}\end{array}}\right.{\omega_{2}}=\left\\{{\begin{array}[]{*{20}{c}}{{{\left({m!}\right)}^{\frac{{-1}}{m}}},}&{{\rm{if\,\,}}m>1}\\\
{1,}&{{\rm{if\,\,}}m<1}\end{array}}\right.,$ (57)

and step $(c)$ is given by the binomial theorem, and it is necessary to assume
that ${m_{{\rm{LoS}}}}$ is an integer.

## References

  * [1] S. Sekander, H. Tabassum, and E. Hossain, “Multi-tier drone architecture for 5G/B5G cellular networks: Challenges, trends, and prospects,” _IEEE Communications Magazine_ , vol. 56, no. 3, pp. 96–103, 2018.
  * [2] B. Li, Z. Fei, and Y. Zhang, “UAV communications for 5G and beyond: Recent advances and future trends,” _IEEE Internet of Things Journal_ , vol. 6, no. 2, pp. 2241–2263, 2018.
  * [3] Z. Lou, A. Elzanaty, and M.-S. Alouini, “Green tethered UAVs for EMF-aware cellular networks,” _IEEE Transactions on Green Communications and Networking_ , vol. 5, no. 4, pp. 1697–1711, 2021.
  * [4] Y. Zeng, R. Zhang, and T. J. Lim, “Wireless communications with unmanned aerial vehicles: Opportunities and challenges,” _IEEE Communications Magazine_ , vol. 54, no. 5, pp. 36–42, 2016.
  * [5] M. Alzenad, M. Z. Shakir, H. Yanikomeroglu, and M.-S. Alouini, “FSO-based vertical backhaul/fronthaul framework for 5G+ wireless networks,” _IEEE Communications Magazine_ , vol. 56, no. 1, pp. 218–224, 2018.
  * [6] M. Kishk, A. Bader, and M.-S. Alouini, “Aerial base station deployment in 6G cellular networks using tethered drones: The mobility and endurance tradeoff,” _IEEE Vehicular Technology Magazine_ , vol. 15, no. 4, pp. 103–111, 2020.
  * [7] Y. Huo, X. Dong, T. Lu, W. Xu, and M. Yuen, “Distributed and multilayer UAV networks for next-generation wireless communication and power transfer: A feasibility study,” _IEEE Internet of Things Journal_ , vol. 6, no. 4, pp. 7103–7115, 2019.
  * [8] D. Wang, S. A. Al-Ahmed, and M. Z. Shakir, “Optimized link distribution schemes for ultrareliable and low-latent communications in multilayer airborne networks,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 9, pp. 5866–5873, 2019.
  * [9] M. Alzenad and H. Yanikomeroglu, “Coverage and rate analysis for vertical heterogeneous networks (VHetNets),” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 12, pp. 5643–5657, 2019.
  * [10] D. Herculea, C. S. Chen, M. Haddad, and V. Capdevielle, “Straight: Stochastic geometry and user history based mobility estimation,” in _Proceedings of the 8th ACM International Workshop on Hot Topics in Planet-scale mobile computing and online social networking_ , 2016, pp. 1–6.
  * [11] M. Haenggi, _Stochastic Geometry for Wireless Networks_. Cambridge University Press, 2012.
  * [12] E. Steinmetz, M. Wildemeersch, T. Q. Quek, and H. Wymeersch, “A stochastic geometry model for vehicular communication near intersections,” in _Globecom Workshops (GC Wkshps)_. IEEE, 2015, pp. 1–6.
  * [13] R. Wang, M. A. Kishk, and M.-S. Alouini, “Stochastic geometry-based low latency routing in massive LEO satellite networks,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 58, no. 5, pp. 3881–3894, 2022.
  * [14] R. Wang, A. Talgat, M. A. Kishk, and M.-S. Alouini, “Conditional contact angle distribution in leo satellite-relayed transmission,” _IEEE Communications Letters_ , pp. 1–1, 2022.
  * [15] N. Okati, T. Riihonen, D. Korpi, I. Angervuori, and R. Wichman, “Downlink coverage and rate analysis of low Earth orbit satellite constellations using stochastic geometry,” _IEEE Transactions on Communications_ , vol. 68, no. 8, pp. 5120–5134, 2020.
  * [16] R. Wang, M. A. Kishk, and M.-S. Alouini, “Evaluating the accuracy of stochastic geometry based models for LEO satellite networks analysis,” _IEEE Communications Letters_ , vol. 26, no. 10, pp. 2440–2444, 2022.
  * [17] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” _IEEE Wireless Communications Letters_ , vol. 3, no. 6, pp. 569–572, 2014.
  * [18] A. Al-Hourani, S. Kandeepan, and A. Jamalipour, “Modeling air-to-ground path loss for low altitude platforms in urban environments,” in _global communications conference_. IEEE, 2014, pp. 2898–2904.
  * [19] X. Wang, H. Zhang, Y. Tian, and V. C. Leung, “Modeling and analysis of aerial base station-assisted cellular networks in finite areas under LoS and NLoS propagation,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 10, pp. 6985–7000, 2018.
  * [20] B. Galkin, J. Kibilda, and L. A. DaSilva, “A stochastic model for UAV networks positioned above demand hotspots in urban environments,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 7, pp. 6985–6996, 2019.
  * [21] W. Bao, Y. Li, and B. Vucetic, “User mobility analysis in disjoint-clustered cooperative wireless networks,” in _Proceedings of the Eighteenth ACM International Symposium on Mobile Ad Hoc Networking and Computing, Los Angeles, USA_ , 2018, pp. 211–220.
  * [22] D. Kim, J. Lee, and T. Q. Quek, “Multi-layer unmanned aerial vehicle networks: Modeling and performance analysis,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 1, pp. 325–339, 2019.
  * [23] W. Bao and B. Liang, “Handoff rate analysis in heterogeneous wireless networks with Poisson and Poisson cluster patterns,” in _Proceedings of the 16th ACM International Symposium on Mobile Ad Hoc Networking and Computing, Paderborn, Germany_ , 2015, pp. 77–86.
  * [24] A. Talgat, M. A. Kishk, and M.-S. Alouini, “Nearest neighbor and contact distance distribution for binomial point process on spherical surfaces,” _IEEE Communications Letters_ , vol. 24, no. 12, pp. 2659–2663, 2020.
  * [25] ——, “Stochastic geometry-based analysis of LEO satellite communication systems,” _IEEE Communications Letters_ , vol. 25, no. 8, pp. 2458–2462, 2021.
  * [26] M. Alzenad, A. El-Keyi, F. Lagum, and H. Yanikomeroglu, “3-D placement of an unmanned aerial vehicle base station (UAV-BS) for energy-efficient maximal coverage,” _IEEE Wireless Communications Letters_ , vol. 6, no. 4, pp. 434–437, 2017.
  * [27] D. Wackerly, W. Mendenhall, and R. L. Scheaffer, _Mathematical Statistics with Applications_. Cengage Learning, 2014\.
  * [28] H. Alzer, “On some inequalities for the incomplete gamma function,” _Mathematics of Computation_ , vol. 66, no. 218, pp. 771–778, 1997.
  * [29] M. N. Kulkarni, S. Singh, and J. G. Andrews, “Coverage and rate trends in dense urban mmWave cellular networks,” in _IEEE Global Communications Conference_. IEEE, 2014, pp. 3809–3814.
  * [30] J. G. Andrews, F. Baccelli, and R. K. Ganti, “A tractable approach to coverage and rate in cellular networks,” _IEEE Transactions on communications_ , vol. 59, no. 11, pp. 3122–3134, 2011.
  * [31] W. U. Mondal, A. A. Sardar, and G. Das, “Economic analysis of cognitive underlay networks: A nash bargaining based approach,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 2, pp. 2024–2029, 2021.
  * [32] R. Wang, M. A. Kishk, and M.-S. Alouini, “Ultra-dense LEO satellite-based communication systems: A novel modeling technique,” _IEEE Communications Magazine_ , vol. 60, no. 4, pp. 25–31, 2022.
  * [33] J. Peng, W. Tang, and H. Zhang, “Directional antennas modelling and coverage analysis of UAV-assisted networks,” _IEEE Wireless Communications Letters_ , 2022.
  * [34] H. Huang and A. V. Savkin, “Navigating UAVs for optimal monitoring of groups of moving pedestrians or vehicles,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 4, pp. 3891–3896, 2021.
  * [35] J. G. Andrews, A. K. Gupta, and H. S. Dhillon, “A primer on cellular network analysis using stochastic geometry,” available online: https://arxiv.org/abs/1604.03183.
  * [36] H. ElSawy, A. Sultan-Salem, M.-S. Alouini, and M. Z. Win, “Modeling and analysis of cellular networks using stochastic geometry: A tutorial,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 1, pp. 167–203, 2016.
  * [37] T. Bai and R. W. Heath, “Coverage and rate analysis for millimeter-wave cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 2, pp. 1100–1114, 2014.

  *[PPP]: Poisson point process
  *[LEO]: low earth orbit
  *[UAVs]: Unmanned Aerial Vehicle
  *[UAV]: Unmanned Aerial Vehicle
  *[PDF]: probability density functions
  *[SINR]: signal-to-interference plus noise ratio
  *[AWGN]: additive white Gaussian noise
  *[CDF]: Cumulative Distribution Function
  *[CCDF]: complementary cumulative distribution function