PSO

     particles swarm optimization
RF

    radio frequency
THz

    terahertz
HRTRN

    hybrid RF and THz relay network
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

#  Terrain-Based UAV Deployment:
Providing Coverage for Outdoor Users

Zhengying Lou, Ruibo Wang, Baha Eddine Youcef Belmekki, Mustafa A. Kishk,
Member, IEEE, and Mohamed-Slim Alouini, Fellow, IEEE Zhengying Lou, Ruibo
Wang, Baha Eddine Youcef Belmekki 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, Maynooth University, Maynooth, W23 F2H6, Ireland.
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>mustafa.kishk @mu.ie;
slim.alouini@kaust.edu.sa).

###### Abstract

Deploying unmanned aerial vehicle (UAV) networks to provide coverage for
outdoor users has attracted great attention during the last decade. However,
outdoor coverage is challenging due to the high mobility of crowds and the
diverse terrain configurations causing building blockage. Most studies use
stochastic channel models to characterize the impact of building blockage on
user performance and do not take into account terrain information. On the
other hand, real-time search methods use terrain information, but they are
only practical when a single UAV serves a single user. In this paper, we put
forward two methods to avoid building blockage in a multi-user system by
collecting prior terrain information and using real-time search. We proposed
four algorithms related to the combinations of the above methods and their
performances are evaluated and compared in different scenarios. By adjusting
the height of the UAV based on terrain information collected before
networking, the performance is significantly enhanced compared to the one when
no terrain information is available. The algorithm based on real-time search
further improves the coverage performance by avoiding the shadow of buildings.
During the execution of the real-time search algorithm, the search distance is
reduced using the collected terrain information.

###### Index Terms:

UAV deployment, terrain, coverage probability, networking, optimization.

##  I Introduction

Achieving ubiquitous coverage for outdoor regions is one of the critical
issues in the fifth generation (5G) and beyond 5G communication networks [1,
2, 3]. Outdoor coverage includes providing connectivity for outdoor residents,
firefighters during a wildfire, rescue teams in emergency situations, and
soldiers during military operations. In particular, non-terrestrial networks
(NTN) platforms, comprising satellite constellations, high altitude platform
stations (HAPS), and low altitude platform stations (LAPS), play an
indispensable role in achieving seamless outdoor coverage [3, 4]. Compared to
conventional ground-based infrastructure, NTN platforms offer swifter
construction and deployment. [5, 6]. However, outdoor coverage is challenging
due to the high mobility of the crowds, the dominance of non-line-of-sight
(NLoS) propagation paths [7], and the diverse terrain configurations [8]. The
self-adaptive and self-evolving micro-networks that are expected in the sixth
generation (6G) vision will provide coverage for outdoor users [9, 10, 11,
12]. Owing to their mobility, unmanned aerial vehicle (UAV) networks are well-
suited for such a highly flexible network architecture [13]. Furthermore,
UAVs’ lower operating altitude grants the benefit of reduced communication
latency, which is vital for applications requiring real-time data
transmission.

For an energy-limited network, the power received by users will directly
determine the communication quality of service [14, 15, 5]. Moreover, the
common metrics used to assess aerial networks, such as coverage probability,
achievable data rate, and energy efficiency, are directly related to the
users’ received power [4, 16, 17, 18]. In scenarios where the UAV provides
coverage for outdoor users, the user’s received power can be enhanced by
shortening the distance between the UAV and the user and by avoiding building
blockage [19, 20]. Most of the research related to UAV networking attaches
importance to the terrain influence such as the blockage of buildings [21].
Although terrain impact on performance is a major focus in UAV networks, most
studies use basic models to describe the impact of building blockage on the
user’s performance. In the following, we present three challenges that impede
the usage of elaborated and sophisticated models.

The first challenge is that lowering the UAV altitude reduces the distance
between the UAV and the user but increases the probability of NLoS
propagation. A widely applied stochastic channel model [22] provides an
acceptable solution to address this challenge. The second challenge is related
to the irregularity of the terrains that can significantly impact the air-to-
ground (ATG) link quality. In this case, a slight movement of the UAV may lead
to a significant change in channel gain [23]. The real-time search approach,
as a new approach emerging in recent years, is leveraged to solve the fine-
grained UAV deployment problem, which addresses the second challenge[24, 25].
The third challenge pertains to simultaneously managing multiple users’
quality of service [26]. When the UAV attempts to establish a direct path with
one user, the original direct path with another user might be blocked.
Consequently, when the number of users being served increases, the difficulty
of deployment rises dramatically. While previous research has made significant
progress in addressing the first two challenges individually, the unique
challenge of building blockage in multi-user systems has remained unexplored,
and strategies to address all three challenges simultaneously have yet to be
provided. In order to tackle these three challenges simultaneously, we have
introduced novel algorithms that leverage both prior terrain information and
real-time search. The fusion of these elements represents a pioneering
approach, setting our work apart from existing research. Importantly, our
approach also incorporates coverage probability estimation and user
classification, further underscoring its innovative nature.

###  I-A Related Works

Most of the existing works are based on a stochastic channel model to estimate
the ATG wireless channel model [27, 28, 29]. The stochastic channel model
answers the following question: how the relative positions between UAVs and
the served users affect the probability of establishing line-of-sight (LoS)
links between them?[22, 30]. In [27], the authors divide ATG propagation modes
into two main types, i.e., LoS and NLoS, while their probability distribution
is modeled as a function of the elevation angle from user to UAV, however,
shadow fading is not considered. Subsequently, the authors [28] provide a
general ATG path loss prediction model based on the environment statistical
parameters related to the receiver’s location. In addition, the authors [22]
propose a simple modified sigmoid function with two environment parameters to
approximate the probability of LoS. Although the statistical model of the ATG
channel can provide an estimation for the expected service level [31], it
cannot exploit the terrain information for UAV positioning.

Compared to UAV deployment based on stochastic channel models, real-time
search methods are more granular and utilize terrain information more
efficiently [24, 25]. In [32], the authors maximize the coverage area of UAVs
while users are located on uneven terrain, i.e., the altitudes of users are
different, but the environment is still coarse-grained. More recently, some
literature combine actual and fine-grained terrain with a UAV-assisted system
in which a UAV acts as a relay for a single user and base station (BS) to
extend the coverage of BS and boost the capacity of a shadowed user [33, 34].
The authors [33] leverage local terrain information to optimize the placement
and trajectory of the UAV. In [34], the authors seek the optimal UAV placement
and demonstrate the superior performance of their proposed terrain-based
approach in a real-world urban topology.

In summary, stochastic channel model-based methods are generally inaccurate
since they use the basic characteristics of the terrain. Although the real-
time search method is more effective in using terrain information, methods in
existing literature only apply to the scenario where one UAV serves one user.
Fortunately, these stochastic channel model-based methods and real-time search
methods are used before and during networking, respectively. As a result, we
are able to combine the benefits of the two approaches, enhance the accuracy
of the former, and expand the use of the latter to accommodate numerous users.

###  I-B Contribution

This paper designs a terrain-based UAV deployment system, divided into four
cases according to the availability of $(I_{1})$ prior terrain information and
$(I_{2})$ real-time search. $(I_{1})$ requires UAV to have sufficient
resources for terrain data collection before providing coverage to users,
while $(I_{2})$ demands UAV to have sufficient electric power and a fast speed
compared to the users for real-time search. The contributions related to each
case are as follows:

  * •

Both $(I_{1})$ and $(I_{2})$ are not available: An improvement of the weighted
$K$-nearest neighbor (WKNN) algorithm, called barycenter-inspired algorithm
(BIA), is proposed. The BIA considers the coverage characteristics of the UAV
network in weight (mass density function) design.

  * •

$(I_{1})$ is available, and $(I_{2})$ is not available: We extract from the
prior terrain information as two terrain-related feature parameters using
ridge regression. The stochastic channel-based positioning algorithm (SCPA) is
proposed based on the parameters for UAV position selection.

  * •

$(I_{1})$ is not available, and $(I_{2})$ is available: We propose the multi-
user-based real-time search algorithm (MRSA), where the UAV can obtain terrain
information and determine the next position while providing network services.
The MRSA applies to the scenario where one UAV serves multiple users with a
linear search length. In particular, the UAV deployment location obtained by
MRSA is proved to be conditionally optimal when a UAV serves two users.

  * •

Both $(I_{1})$ and $(I_{2})$ are available: With prior terrain information, we
propose a resident classification method to narrow the range of users that UAV
needs to serve. Based on this method, we improve MRSA and further propose the
hybrid deployment algorithm (HDA). Compared to MRSA, HDA requires a shorter
search distance and can find the UAV deployment location with better coverage
performance.

Please note that MRSA and HDA are algorithms employed for the first time in
multi-user UAV systems. Additionally, what distinguishes them from existing
real-time search algorithms is the integration of probability estimation and
user classification.

Figure 1: Illustration of the system model configurations.

##  II System Model

In this section, we provide a scenario description, channel model, and
optimization objectives related to the study.

###  II-A Scenario Description

We consider a suburban region containing several buildings, one UAV, and $K$
outdoor ground users. Users are required to upload their locations when
requesting network access services. The UAV dynamically adjusts its position
based on the users’ locations. To avoid collisions with buildings, the minimum
height of the UAV, denoted by $h_{\min}$, is higher than the tallest building
as depicted in Fig 1. The above assumptions are common and reasonable in UAV
networking [35]. However, the specific details of the terrain are unknown to
the UAV since most of the available databases do not provide detailed terrain
information, such as the height and shape of buildings as well as other signal
barriers such as large vehicles. Finally, we assume that there are no other
public access points (APs) in the area, and outdoor users can only access the
network via the UAV. Then, the UAV establishes a dependable LoS connection
with the core network to meet the backhaul requirements, owing to the higher
altitudes at which both the UAV and base station are typically deployed
compared to the surrounding buildings.

###  II-B Channel Model

The signal received by the user is considered in a LoS when no building is
blocking the signal between the user and the UAV. Otherwise, the signal is
considered in a NLoS. The UAV-user channel undergoes large-scale and small-
scale fading. Therefore, the received power $S_{Q}$,
$Q=\left\\{\rm{LoS},\rm{NLoS}\right\\}$ is expressed as

$S_{Q}\left(r\right)=\eta_{Q}\,\zeta\,G_{Q}r^{-\alpha_{Q}},$ (1)

where $\zeta$ is the transmit power, $\eta_{Q}$ is the mean additional losses
[36], and $r$ is the Euclidean distance between the user and the UAV. The
large-scale fading is denoted by $\eta_{Q}\,r^{-\alpha_{Q}}$, where
$\alpha_{Q}$ is the path-loss exponent. Considering the significant multipath
effect of buildings, the small-scale fading $G_{Q}$ is modeled by a
Nakagami-$m$ fading channel with a shape parameter and scale parameter
$(m_{Q},\frac{1}{m_{Q}})$. The probability density functions (PDF) of the
power gains $G_{Q}$ is given by [37]

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

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

Therefore, the overall coverage probability, which is the probability that the
received signal-to-noise ratio (SNR) is greater than a predefined threshold
$\gamma$, can be derived by

$P_{C}\left(h;r_{1},r_{2},...,r_{K}\right)=\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}\left[\frac{S_{Q,k}\left(r_{k}\right)}{\sigma^{2}}>\gamma\right],$
(3)

where $h$ is the height of the UAV, $r_{k}$ is the Euclidean distance between
the $k^{th}$ user and the UAV, $\sigma^{2}$ is the environmental noise
received by the user, and $Q$ is either LoS or NLoS for a given user $k$.

###  II-C Optimization Objectives

Our goal is to determine the UAV position that maximizes the overall coverage
probability $P_{C}$, which we define as the optimal position.

######  Definition 1 (Terrain-based Optimal Position).

The position that has the maximum overall coverage probability within the
three-dimensional region accessible to the UAV is called a terrain-based
optimal position (simply denoted as optimal).

Considering the above channel model, the optimal position can be determined
using the following two approaches: reducing the distance $r_{k}$, and
establishing LoS link rather than NLoS link with user $k$. Then, our
optimization problem can be formulated as

$\displaystyle\underset{\boldsymbol{h;r_{k},1\leq k\leq K}}{\text{maximize}}$
$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}\left[\frac{S_{Q,k}\left(r_{k}\right)}{\sigma^{2}}>\gamma\right],$
(4a) subject to: $\displaystyle h>h_{\min},$ (4b)

where the objective function is the overall coverage probability that can be
seen as a weighted sum of users’ coverage probability such that the weight of
each user is the same. The height of the UAV $h$, and the Euclidean distance
between the UAV and each user $r_{k},1\leq k\leq K$, are treated as
optimization variables. The only constraint (4b) is enforced to prevent
collisions between the UAV and buildings. Finally, it is noted that the
position with the highest instant coverage probability or received signal
power is not necessarily the optimal position on the trajectory because of the
small-scale fading fluctuations. Therefore, the UAV can only utilize the
instant received signal power to check the blockage in real-time search.

##  III The Acquisition and Processing of Terrain Information

###  III-A LoS Probability Estimation

In this subsection, we first introduce the LoS probability expression with
undetermined parameters, then proceed to describe how these parameters can be
estimated according to the data collected by UAV. Finally, we discuss the LoS
probability in more detail.

We consider that the UAV can obtain the characteristics of the LoS probability
by collecting a small amount of data in advance, but such data is not enough
for reconstructing the terrain. The authors in [22, 31] express the
probability of establishing a LoS link as a modified sigmoid function

$P_{\rm{LoS}}\left(h,r\right)=\frac{1}{{1+a\exp\left({-b\left({\frac{{180}}{\pi}{{\tan}^{-1}}\left({\frac{{{h}}}{\sqrt{r^{2}-h^{2}}}}\right){-a}}\right)}\right)}}.$
(5)

Based on the above function, the relationship between
$P_{\rm{LoS}}\left(h,r\right)$ and $h,r$ can be determined by two terrain-
related parameters, $a$ and $b$. The values of $a$ and $b$ are related to the
density, height, and floor area of the buildings. Since the establishment of
LoS link and NLoS link are mutually exclusive events, the corresponding two
probabilities are complementary, that is,
$P_{{\rm{NLoS}}}\left(h,r\right)=1-P_{{\rm{LoS}}}\left(h,r\right).$

Next, we explain how to collect the data set and estimate parameters $a$ and
$b$ based on the data set. In order to facilitate data collection, we define
the elevation angle $\theta$ as

$\theta=\frac{{180}}{\pi}{{\tan}^{-1}}\left({\frac{{{h}}}{\sqrt{r^{2}-h^{2}}}}\right).$
(6)

For the data collection of a fixed $\theta=\theta_{i}$, the height of the UAV
$h$ is adjusted to satisfy the relationship between $h$, $\theta_{i}$ and the
horizontal distance from the UAV to a randomly chosen user
$\sqrt{r^{2}-h^{2}}$ in (6). In order to ensure fairness during the data
collection process, the UAV modifies its position and associates various users
while maintaining a fixed $\theta_{i}$. Then, the ratio between the number of
samples of LoS link established and the total number of samples is obtained,
denoted as $t_{i}$. Finally, the above steps are repeated by changing $\theta$
to obtain a data set containing $N$ tuples:
$\mathcal{T}=\\{(t_{1},\theta_{1}),(t_{2},\theta_{2}),\dots,(t_{N},\theta_{N})\\}$.

Based on the data set $\mathcal{T}$, we choose the mean square error (MSE)
function as the cost function. Regularization is added in order to enhance the
generalization ability and to reduce the influence of the edge effect (non-
uniform sampling in edge regions due to limited map size). Therefore, the
parameter estimation problem is described by the following nonlinear ridge
regression problem

$\left(a^{*},b^{*}\right)=\underset{{a,b}}{\arg\min}\sum_{i=1}^{N}\left(t_{i}-\frac{1}{{1+a\exp\left({-b\left({\theta_{i}-a}\right)}\right)}}\right)^{2}+\lambda_{1}\left(a-\widehat{a}\right)^{2}+\lambda_{2}\left(b-\widehat{b}\right)^{2},$
(7)

where $\widehat{a}$ and $\widehat{b}$ are empirical parameters that provided
in [38], $\lambda_{1}$ and $\lambda_{2}$ are regularization coefficients.
Since $\widehat{a}$ is much larger than $\widehat{b}$, $\lambda_{2}$ should
also be much larger than $\lambda_{1}$ to ensure that both arguments have
similar learning rates.

Finally, we add two relevant comments to the LoS probability estimation
problem. Firstly, the modified sigmoid function in (5) may not be the only
choice, and good estimates can be made by improving other commonly used
machine learning functions. We improve the Tanh and ReLu functions and also
use only two parameters to estimate the LoS probability. Secondly, we apply
the trust region reflective algorithm to solve the optimization problem in
(7), and the initial point is set as $(\widehat{a},\widehat{b})$.

###  III-B Coverage Probability and Resident Classification

In this subsection, we first give the motivation for classification. Then, a
closed-form expression for the coverage probability is derived. Finally, we
show how the users are classified based on the coverage probability.

The authors in [39] introduce a probabilistic node sensing model and divide
the receivers into three classes according to their distances from the
transmitter: definitely covered ($C_{1}$), probably covered ($C_{2}$), and
cannot be covered ($C_{3}$). Through the classification, the UAV should focus
on $C_{2}$ users to improve the overall coverage, because it does not need to
focus on $C_{1}$ (close) users and the $C_{3}$ (distance) users, since the
former are already covered, the latter cannot be covered.

The second motivation for user classification is that the coverage probability
of $C_{2}$ users decreases almost linearly with the increase of the Euclidean
distance between the UAV and the user. Using this important property, the
expression of overall coverage probability can be simplified in some cases.
Based on (1) and (5), the coverage probability is given in the following
theorem.

######  Lemma 1.

Given that the height of UAV is $h$ and distance between the $k^{th}$ user and
the UAV is $r_{k}$, the coverage probability of the $k^{th}$ user is given by

$P_{C,k}\left(h,r_{k}\right)=P_{\rm{LoS}}\left(h,r_{k}\right)P_{C,k}^{\rm{LoS}}\left(r_{k}\right)+P_{\rm{NLoS}}P_{C,k}^{\rm{NLoS}}\left(r_{k}\right),$
(8)

where $P_{\rm{LoS}}\left(h,r_{k}\right)$ is defined in (5), and
$P_{{\rm{NLoS}}}\left(h,r_{k}\right)=1-P_{{\rm{LoS}}}\left(h,r_{k}\right).$
$P_{C,k}^{Q}\left(h,r_{k}\right)$ is the coverage probability under the
condition that $Q$ link is established, where
$Q=\left\\{\rm{LoS},\rm{NLoS}\right\\}$. The expression of
$P_{C,k}^{Q}\left(h,r_{k}\right)$ is given by

$\begin{split}P_{C,k}^{Q}\left(r_{k}\right)=\exp\left(-\mu_{Q}\left(r_{k}\right)\right)\sum_{n=0}^{m_{Q}-1}\frac{\mu_{Q}\left(r_{k}\right)^{n}}{n!},\end{split}$
(9)

where $m_{Q}$ is shape parameter of Nakagami-$m$ fading and
$\mu_{Q}\left(r_{k}\right)$ is given by

$\mu_{Q}\left(r_{k}\right)=\eta_{Q}^{-1}\zeta^{-1}\gamma\,r_{k}^{\alpha_{Q}}\sigma^{2}.$
(10)

###### Proof.

See Appendix A. ∎

The conditional coverage probability expression of in (9) contains the
exponential term
$\exp\left(-\mu_{Q}\left(r_{k}\right)\right)$, therefore, when
$\mu_{Q}\left(r_{k}\right)\gg 1$, the conditional coverage probability
converges rapidly to 0; whereas when $0<\mu_{Q}\left(r_{k}\right)\ll 1$, the
conditional coverage probability converges rapidely to 1. This clarifies the
reasoning behind the above classification. Only a small range of $r_{k}$ makes
$\mu_{Q}\left(r_{k}\right)\approx 1$ and brings a significant change in the
conditional coverage probability.

Note that the coverage probability is related to the LoS probability obtained
in subsection III-A, while the conditional coverage probability is not.
Accordingly, we denote by $\varepsilon$ the degree of classifications, and we
provide the following two types of classifications:

  * •

Non-terrain-based classification: When
$P_{C,k}^{\rm{NLoS}}\left(r_{k}\right)>1-\varepsilon$, user $k$ belongs to
$C_{1}$; when $P_{C,k}^{\rm{LoS}}\left(r_{k}\right)<\varepsilon$, user $k$
belongs to class $C_{3}$; otherwise user $k$ is classified into the class
$C_{2}$.

  * •

Terrain-based classification: When $P_{C,k}\left(r_{k}\right)>1-\varepsilon$,
user $k$ belongs to $C_{1}$; when $P_{C,k}\left(r_{k}\right)<\varepsilon$,
user $k$ belongs to class $C_{3}$; otherwise when $\varepsilon\leq
P_{C,k}\left(r_{k}\right)\leq 1-\varepsilon$, user $k$ is classified into the
class $C_{2}$.

###  III-C Real-Time Search in Simple Scenarios

In this subsection, we will start with the real-time terrain-based search for
two simple scenarios. The conclusion obtained from the following analysis of
simple scenarios can inspire the design of the MRSA and HDA.

  * •

Scenario I: There are no buildings, and the height of the UAV’s optimal
position is, therefore, the lowest altitude $h_{\min}$. This scenario is
analyzed in order to find a feasible optimization objective.

  * •

Scenario II: The UAV provides coverage for only two users. Based on the
analysis of scenario II, a practical search strategy is proposed.

For the above two scenarios, the UAV searches along a trajectory specified by
a strategy $\mathscr{L}$ and chooses a point on the trajectory as the final
positioning result. For the above search methods (MRSA and HDA), the two
characteristics that we are most interested in are optimality, defined in
section II, and linearity.

######  Definition 2 (Linearity of the Trajectory).

Linearity is a function of (i) the distance among the users, and (ii) the UAV
height. For example, the distance among the users increases, (the upper bound
of) the length of the trajectory of $\mathscr{L}$ increases linearly, in that
case, the strategy $\mathscr{L}$ is called a linear strategy.

Furthermore, from both latency and energy perspectives, the length of
trajectory is highly limited. Therefore, the linear search length of the UAV
is necessary. Compared to the flight duration of the UAV, the computational
delay is significantly shorter, so we place a greater emphasis on the analysis
of UAV search length. Combining linearity and optimality, the following
definition is given.

######  Definition 3 (Linear Optimal Strategy).

If a linear strategy $\mathscr{L}$ can always pass the optimal position, then
we call $\mathscr{L}$ a linear optimal strategy.

A squared optimal strategy, which involves searching the whole plane to find
the optimal solution, exists in scenario I. However, whether a linear optimal
strategy, which aims to find the optimal solution by searching the space
incrementally with a linear growth pattern, exists in scenario I is discussed
in the following proposition.

######  Proposition 1.

In scenario I, the relationship between the existence of linear optimal
trajectory $\mathscr{L}$ and the number of users is as follows:

  * •

When the UAV provides coverage to one or two users, a linear optimal
trajectory always exists.

  * •

When the UAV provides coverage to three users, an approximate linear optimal
trajectory can be found by approximating the problem of searching the optimal
position to the problem of finding the Fermat-Weber point [40].

  * •

When the UAV provides coverage to more than four users, the linear optimal
trajectory does not exist.

###### Proof.

See Appendix B. ∎

Although the above analysis can provide inspiration for the real-time terrain-
based search algorithm, it is challenging to provide effective search
trajectories for large numbers of users. Therefore, we define the average SNR
as

$\overline{\rm{SNR}}_{Q}(r)=\mathbb{E}_{G}\left[\frac{S_{Q}(r)}{\sigma^{2}}\right]=\frac{\eta_{Q}\zeta
r^{\alpha_{Q}}}{\sigma^{2}},$ (11)

where $Q=\\{\rm{LoS},\rm{NLoS}\\}$. Note that optimizing other performance
metrics, such as average data rate and energy efficiency [41], ultimately
boils down to optimizing the average SNR. While our focus lies on optimizing
overall coverage probability, it indirectly impacts and improves these
associated metrics. The following definition is given to simplify the
optimization objective by seeking the linear suboptimal trajectory instead of
the linear optimal one.

######  Definition 4 (Linear $\gamma$-suboptimal Strategy).

A strategy $\mathscr{L}$ is said to be linear $\gamma$-suboptimal (simply
denoted as $\gamma$-suboptimal) if (i) there is one and only one
$\gamma$-suboptimal position that satisfies the condition that the average SNR
from this position to any user in $C_{2}$ is greater than or equal to
$\gamma$, (ii) $\mathscr{L}$ can find the $\gamma$-suboptimal position, and
(iii) $\mathscr{L}$ is a linear strategy.

######  Proposition 2.

In scenario I, there always exists a $\gamma$-suboptimal strategy
$\mathscr{L}$.

###### Proof.

See Appendix C. ∎

Next, we will propose a $\gamma$-suboptimal strategy in scenario II. Since the
UAV provides coverage for only two users, we establish a cylindrical
coordinate system. Take the midpoint of the two users as the origin, the line
passing through the two users as $z$-axis, and the direction perpendicular to
the ground as the direction of $\theta=\frac{\pi}{2}$. Therefore, the users’
cylindrical coordinates are $(\rho,\theta,z)=(0,0,\pm\frac{d}{2})$, where $d$
is the distance between two users.

Figure 2: Heat maps of smaller average SNR between two users from different
perspectives.

Fig. 2 shows heat maps of smaller average SNR between two users from different
perspectives. There are some distinct boundaries in the figure, indicating
that one user switches from not being blocked to being blocked, while the
other user is never blocked by the building. From Fig. 2, we know that only in
the $z=0$ plane, the LoS region is continuous because any point below a point
in the NLoS region also belongs to the NLoS region. Given the special property
of $z=0$ plane, the following strategy is proposed: (i) search on the $z=0$
plane, (ii) decrease $\rho$ until at least one user’s link to the UAV is
blocked by the building and (iii) change the value of $\theta$ until both
links (between the UAV and the two users) are in LoS.

Algorithm 1 Real-Time Search Algorithm for Two Users

1: Input: The starting coordinates of the UAV $(\rho_{u},\theta_{u},0)$,
search granularity $\delta$.

2: Initialize: Boolean variable $\mathcal{B}\leftarrow 0$ if both of the users
establish LoS links with the UAV, otherwise $\mathcal{B}\leftarrow 1$. Let
$z_{u}\leftarrow 0$.

3: while $\mathcal{B}=1$ do

4: $\rho_{u}\leftarrow\rho_{u}+\delta$, and renew the value of $\mathcal{B}$.

5: end while

6: while $\rho_{u}\cos\theta_{u}>h_{\min}$ do

7: if $\mathcal{B}=1$ then

8: $\rho^{*}\leftarrow\rho_{u}$, $\theta^{*}\leftarrow\theta_{u}$,
$\rho_{u}\leftarrow\rho_{u}-\delta$.

9: else

10: $\theta_{u}\leftarrow\theta_{u}-2\arcsin\frac{\delta}{2\rho_{u}}$.

11: end if

12: Renew the value of $\mathcal{B}$.

13: end while

14: $\rho_{u}\leftarrow\rho^{*}$, $\theta_{u}\leftarrow\theta^{*}$.

15: while $\rho_{u}\cos\theta_{u}>h_{\min}$ do

16: if $\mathcal{B}=1$ then

17: $\rho^{*}\leftarrow\rho_{u}$, $\theta^{*}\leftarrow\theta_{u}$,
$\rho_{u}\leftarrow\rho_{u}-\delta$.

18: else

19: $\theta_{u}\leftarrow\theta_{u}+2\arcsin\frac{\delta}{2\rho_{u}}$.

20: end if

21: Renew the value of $\mathcal{B}$.

22: end while

23: if
$\overline{\rm{SNR}}_{\rm{NLoS}}\left(\sqrt{h_{\min}^{2}+\frac{d^{2}}{4}}\right)>\overline{\rm{SNR}}_{\rm{LoS}}\left(\sqrt{{\rho^{*}}^{2}+\frac{d^{2}}{4}}\right)$
then

24:
$\gamma\leftarrow\overline{\rm{SNR}}_{\rm{NLoS}}\left(\sqrt{h_{\min}^{2}+\frac{d^{2}}{4}}\right)$,
$\rho^{*}\leftarrow h_{\min}$, $\theta^{*}\leftarrow 0$.

25: else

26:
$\gamma\leftarrow\overline{\rm{SNR}}_{\rm{LoS}}\left(\sqrt{{\rho^{*}}^{2}+\frac{d^{2}}{4}}\right)$.

27: end if

28: Output: Position ($\rho^{*},\theta^{*},0$) and $\gamma$.

Steps (3) - (6) of Algorithm 1 ensure that the starting position is in the LoS
region for both users. Steps (7) - (15) search the left branch, and (16) -
(25) search the right branch.  The steps (6) and (15) respectively verify
whether the UAV is below $h_{\min}$ in the left branch search and the right
branch search.

As for the final position of the UAV search trajectory, there are two
potential $\gamma$-suboptimal positions in the LoS and NLoS regions. Steps
(23) - (27) compare these two positions. The potential $\gamma$-suboptimal
position in the NLoS region is $(h_{\min},0,0)$ (given that the position
$(h_{\min},0,0)$ is not located in the common LoS region of both users). In
the LoS region, the potential $\gamma$-suboptimal position is the last
position in the trajectory that the UAV stays in the common LoS region of both
users.

Even though the average SNR is not available to the UAV, it can infer whether
the users are blocked or not.  Furthermore, for high-rise building scenarios
[28], there may be situations where the UAV’s altitude becomes too high. In
such cases, we can set an upper limit for $\rho_{u}$. During the execution of
step (4), if this limit is exceeded, the algorithm terminates and outputs
$(\frac{d_{\min}}{\cos\theta_{u}},\theta_{u},0)$ as the deployment position
for the UAV.

######  Proposition 3.

In scenario II, we have the following statements for Algorithm 1 and its
output $(\rho^{*},\theta^{*},0)$:

  * •

When $2\rho^{*}\leq d$, the strategy proposed by Algorithm 1 is
$\gamma$-suboptimal. In other words, $(\rho^{*},\theta^{*},0)$ is the only
position where the average SNRs of both users are greater than or equal to
$\gamma$.

  * •

When $2\rho^{*}>d$, the strategy proposed by Algorithm 1 is not
$\gamma$-suboptimal and
$\gamma<\overline{\rm{SNR}}_{\rm{LoS}}\left(\frac{2\rho^{*}d}{\sqrt{\left(2\rho^{*}\right)^{2}+d^{2}}}\right)$.

###### Proof.

See Appendix D. ∎

##  IV Typical Scenarios Algorithms Design

In this section, four typical scenarios and the corresponding four positioning
algorithms are shown in Table I. Prior information available indicates that
the UAV is allowed to collect data and determine the prior probability of
establishing LoS link in advance. Since the users move slower than UAV, the
latter has a time span to search for the best location based on the terrain.
Specifically, at each time interval, based on users’ locations, the UAV
adjusts its own location by recalculating its position using the corresponding
algorithm in different scenarios. At last, the granularity $\delta$ is defined
for the positioning accuracy.

TABLE I: Positioning algorithms for four typical scenarios. Algorithm |  | Prior information required
---
| Terrain information utilized
---
| Real-time search required
---
BIA | No | No | No
| SCPA
---
Yes | Yes | No
MRSA | No | Yes | Yes
HDA | Yes | Yes | Yes

###  IV-A Barycenter-Inspired Algorithm

The Algorithm 2, that is, the BIA is designed for the scenario when prior
information is not available. It can be regarded as a variant of the
clustering algorithm.

We first establish a Cartesian coordinate system with the origin on the ground
and the $k^{th}$ user’s location can be expressed as
$\left(x_{k}^{r},y_{k}^{r},0\right)$. After that, we give different weights
$\Lambda\left(r\right)$ to users who are at different distances from the UAV,
and the weight $\Lambda\left(r\right)$ is called the mass density. In each
iteration, the UAV takes the user’s barycenter as its horizontal position.
After adjusting the position of the UAV, the weight $\Lambda\left(r\right)$
also changes, which further leads to the change of barycenter. Therefore, the
position of the UAV needs to be updated. According to the above deceptions,
the BIA is given as follows.

Algorithm 2 Barycenter-inspired algorithm (BIA).

1: Input: Locations of users $\left(x_{k}^{r},y_{k}^{r},0\right)$, mass
density function $\Lambda\left(r\right)$, maximum number of iterations $N$ and
granularity $\delta$.

2: Initialize: deploy the UAV at
$\left(x_{0}^{u},y_{0}^{u},z_{0}^{u}\right)\leftarrow\left(\sum_{k=1}^{K}\frac{x_{k}^{r}}{K},\sum_{k=1}^{K}\frac{y_{k}^{r}}{K},h\right)$,
initialize $n\leftarrow 1$.

3: repeat

4: for $k=1:K$ do

5:
$r_{k}\leftarrow\sqrt{\left(x_{n-1}^{u}-x_{k}^{r}\right)^{2}+\left(y_{n-1}^{u}-y_{k}^{r}\right)^{2}+h^{2}}$.

6: end for

7:
$\left(x_{n}^{u},y_{n}^{u},z_{n}^{u}\right)\leftarrow\left(\frac{\sum_{k=1}^{K}\Lambda\left(r_{k}\right)x_{k}^{r}}{\sum_{k=1}^{K}\Lambda\left(r_{k}\right)},\frac{\sum_{k=1}^{K}\Lambda\left(r_{k}\right)y_{k}^{r}}{\sum_{k=1}^{K}\Lambda\left(r_{k}\right)},h\right)$.

8: $n\leftarrow n+1$.

9: until $n\geq N\ \textbf{or}{}\
\sqrt{\left(x_{n}^{u}-x_{n-1}^{u}\right)^{2}+\left(y_{n}^{u}-y_{n-1}^{u}\right)^{2}}\leq\delta$.

10: Output: The final location of UAV
$\left(x_{n}^{u},y_{n}^{u},z_{n}^{u}\right)$.

There are three remarks on the BIA. Firstly, since the height and distribution
of buildings are unknown to the UAV, the flight height $h$ cannot be optimized
by the distribution of users. Therefore, the algorithm only deploys the
horizontal position of UAV.  Secondly, BIA and Algorithm 1 are expressed by
different coordinate systems. Since BIA is an improvement of the WKNN
algorithm, it is more suitable to be presented in a Cartesian coordinate
system. In order to facilitate the representation of contour lines in
Proposition 3 and the proof of $\gamma$-suboptimal positions, Algorithm 1
utilizes a cylindrical coordinate system for modeling. Finally, we compare
several segmented mass density functions $\Lambda(r)$ in section V.

###  IV-B Stochastic Channel-Based Positioning Algorithm

The SCPA is designed for the scenario when prior information (the LoS
probability distribution) is available. Given the assumption that the
computing power is sufficient, the SCPA is designed as an improvement of the
brute force method.

The algorithm set the mean value of user coordinates as the initial point.
Next, the algorithm searches for the position with the maximum total coverage
probability $\sum_{k=1}^{K}P_{C,k}\left(h,r_{k}\right)$ around the initial
point through brute force. Due to the continuity of user movement, the
position of the UAV at this moment is usually close to the position of the
previous moment. Therefore, the initial point can also be set to the position
of the previous moment, and the search region of the brute force method does
not need to be large.

The similarity between the above two algorithms is that the SCPA can be
regarded as the BIA with a different and inexplicit nonlinear mass density
function. The difference between the two algorithms is that the SCPA optimizes
the height of UAV.

Figure 3: Flow chart of algorithms required for MRSA and HDA.

###  IV-C Multi-user-based Real-time Search Algorithm and Hybrid Deployment
Algorithm

In this subsection, we first summarize the search strategy of the simple
scenarios described in subsection III-C, and then extend the above strategies
into a general one. Since the linear optimal strategy is difficult to find
even in simple scenarios, we focus on the linear $\gamma$-suboptimal strategy
and provide the following deductions.

  * •

Deduction from Proposition 1: When there are less than three users in $C_{2}$,
the Fermat-Weber point of users with the height of $h_{\min}$ can be chosen as
the potential optimal position, regardless of the building blockage.

  * •

Deduction from Proposition 2 and Appendix C: The minimum enclosing circle can
be determined by first finding the two furthest apart users in $C_{2}$, then
repeatedly shrinking the circle and replacing the boundary users. The center
of the smallest enveloping circle is taken as the possible horizontal position
of the UAV, regardless of the blockage of buildings.

  * •

Deduction from Proposition 3: In most cases, the trajectory proposed by
Algorithm 1 is able to avoid the blockage and find the $\gamma$-suboptimal
position when there are only two users. Even if the $\gamma$-suboptimal
position cannot be found, the algorithm can gradually improve the coverage of
two users during searching and provide a possible upper bound for $\gamma$.

Based on the above analysis, we propose algorithms under real-time search,
which apply to a general scenario, that is, at least three users need
coverage. The algorithms under real-time search contain the MRSA and HDA, and
the differences between them are reflected in green steps shown in Fig. 3. To
obtain the center of classification $\left(x_{0},y_{0},z_{0}\right)$, MRSA and
HDA will execute BIA and SCPA, respectively. In user classification, the MRSA
and HDA will apply non-terrain-based classification and terrain-based
classification, respectively. Then the second part (yellow steps in Fig. 3)
are inspired by Proposition 2, and the third part (blue steps in Fig. 3) are
inspired by Proposition 3. Finally, when only one or two users require
networking, the above algorithm can be substituted by hovering directly above
the users (inspired by Proposition 1) or by executing Algorithm 1.

Finally, since the impact of $\varepsilon$ on the coverage performance of MRSA
and HDA is difficult to express directly through explicit expressions, we
provide the following criteria for determining the values of $\varepsilon$.

  * •

When the steep slope portion of the curve that depicts coverage probability as
a function of the UAV-user distance is not categorized as $C_{2}$ (examples
are given in Fig. 11), the value of $\varepsilon$ should be reduced.

  * •

When the search trajectory becomes lengthy, causing the UAV to be unable to
perform a round of MRSA or HDA before the next update of user positions, the
value of $\varepsilon$ should be increased.

  * •

When a large number of users are classified as $C_{2}$, making it challenging
for the UAV to locate the common LoS region of these users, the value of
$\varepsilon$ should be increased.

##  V Numerical Results

In this section, we analyze the characteristics of the four algorithms and
compare their performances. In each simulation, we run the code for $10^{4}$
rounds. In each round, we first generate the user’s location based on the
Poisson point process (PPP) within an area of $300~{}$m$~{}\times~{}300~{}$m,
then we use one of the four algorithms to obtain the position of the UAV, and
finally, we place the UAV at this position and calculate the proportion of the
user under coverage. In the simulation, the buildings are modeled as
prisms111The latitude and longitude of the buildings’ location can be
downloaded from the website www.openstreetmap.org. This paper also converts
them into coordinates in the Cartesian coordinate system., and their heights
follow a Rayleigh probability distribution with suburban’s parameter [22].
Note that we have simplified the representation of buildings in simulation due
to the intricacy of actual terrain construction, but in the majority of real-
world scenarios, our algorithm remains effective in the presence of inclined
or irregular shapes buildings. For the sake of fairness, the algorithms are
compared based on the same value of granularity $\delta=1~{}$m. For the sake
of convenience, the simulation parameters are provided in Table II, unless
specified otherwise.

TABLE II: System configuration. Parameter | Symbol | Value
---|---|---
Transmit power | $\zeta$ |  $30$ dBm
Environmental noise power | $\sigma^{2}$ |  $-98$ dBm
Threshold for SNR | $\gamma$ |  $22$ dB
Path-loss exponents of LoS and NLoS |  $\alpha_{\rm{LoS}}$, $\alpha_{\rm{NLoS}}$ |  $2$, $2.3$
Nakagami-$m$ fading shape parameters of LoS and NLoS |  $m_{\rm{LoS}}$, $m_{\rm{NLoS}}$ |  $2$, $1$
Mean additional losses of LoS and NLoS |  $\eta_{\rm{LoS}}$, $\eta_{\rm{NLoS}}$ |  $-35$ dB, $-48$ dB

Figure 4: Diagram of mass density functions when $h=20~{}$m.

Figure 5: A sampling of the BIA.

###  V-A Performance of the Barycenter-Inspired Algorithm

We start with the mass density function and the BIA. When $h=20~{}$m, the
diagram of four mass density functions $\Lambda(r)$ of the BIA is provided in
Fig.5. For the selection of uniform mass density, it can be regarded that the
algorithm will stop at the initialization step. Considering that users who are
close to the UAV are likely to be covered, and users who are far away are
difficult to be covered by small adjustments to the UAV’s position, neither of
these users deserves high weight. According to the above analysis, the
expressions of the remaining three mass density functions are designed and
given in Table III. $R_{\min}$ and $R_{\max}$ in Table III are non-terrain-
based boundaries of $C_{1}$-$C_{2}$ and $C_{2}$-$C_{3}$ in subsection III-B.
In our scenario, they can be observed from Fig. 11 that $R_{\min}=40~{}$m and
$R_{\max}=126~{}$m.

TABLE III: Expressions of three mass density functions. |  | $r\leq\max\\{h,R_{\min}\\}$
---
or $R_{\max}<r$
| $\max\\{h,R_{\min}\\}<r$
---
and $r\leq\frac{1}{2}\sqrt{R_{\max}^{2}+3h^{2}}$
$\frac{1}{2}\sqrt{R_{\max}^{2}+3h^{2}}<r\leq R_{\max}$
Ascending trapezoid | $\Lambda(r)=0$ | $\Lambda(r)=\sqrt{r^{2}-h^{2}}$ | $\Lambda(r)=\frac{1}{2}\sqrt{R_{\max}^{2}-h^{2}}$
Descending trapezoid | $\Lambda(r)=0$ | $\Lambda(r)=\frac{1}{2}\sqrt{R_{\max}^{2}-h^{2}}$ | $\Lambda(r)=\sqrt{R_{\max}^{2}-h^{2}}-\sqrt{r^{2}-h^{2}}$
Triangular | $\Lambda(r)=0$ | $\Lambda(r)=\sqrt{r^{2}-h^{2}}$ | $\Lambda(r)=\sqrt{R_{\max}^{2}-h^{2}}-\sqrt{r^{2}-h^{2}}$

Fig. 5 shows a typical sampling of the BIA. Except for the uniform mass
density function, the horizontal positions of the UAV obtained by the
remaining three mass density functions are usually close to each other,
regardless of the height for $h=20~{}$m or $h=60~{}$m.

Figure 6: Influences of UAV height and mass density function on coverage
probability.

Figure 7: Parameter estimation of LoS probability distribution.

Fig. 7 illustrates that deploying the UAV at relatively low heights and
applying non-uniform density functions can effectively improve coverage
performance.  Note that the coverage probability varies at different locations
due to the varying distances to the UAV and the impact of building blockage.
After collecting coverage probability samples from different locations across
the entire area, we can obtain the cumulative distribution function (CDF) of
coverage probability, which is the $y$-axis of Fig. 7 (as well as Fig. 9, 13
and, 15). When the UAV is deployed at 20 m, the coverage probability is
between 10% and 50%. When the UAV is deployed at 60 m, the coverage
probability is almost always less than 40%. By comparing CDFs from different
density functions at the same height, it can be seen that it is more effective
in reducing the weight of users who are farther from the UAV than those who
are closer.

###  V-B Performance of Stochastic Channel-Based Positioning Algorithm

In this subsection, we first provide the numerical results pertaining to the
estimation of the LoS probability distribution. Then, we analyze the
performance of SCBA.

TABLE IV: Comparison of parameters and MSEs under different estimates. | $\lambda_{1}$ | $\lambda_{1}$ | $a$ | $b$ | MSE
---|---|---|---|---|---
Empirical parameters | N/A | N/A | 4.88 | 0.43 | 0.0220
Regularization I | 0.001 | 0.1 | 4.52 | 0.17 | 0.0074
Regularization II | 0.01 | 0.01 | 4.76 | 0.15 | 0.0071
Non-regularization | 0 | 0 | 1.93 | 0.07 | 0.0017

Firstly, the UAV obtains the LoS probability under different elevation angles
by random sampling. The probability is estimated as the proportion of LoS
links established between the user and the UAV, which is shown as the red dots
in Fig. 7. Based on these data, we present the results of regularized and non-
regularized nonlinear regression in Fig. 7. Empirical parameters of sub-urban
environments ($\widehat{a}=4.88$, $\widehat{b}=0.43$, provided by authors in
[38]) are used for comparison and regularization.

Table IV provides the parameters $a$ and $b$ and the mean square error (MSE)
obtained by different estimations. Both the MSEs in Table IV and the curves in
Fig. 7 show that the empirical parameters do not fit well with the LoS
probability distribution. In fact, since only two parameters need to be
estimated, relatively accurate results can be obtained by only a small data
set with hundreds of samples. Note that if dozens of data are collected for
only several elevation angles, smoothing is necessary to eliminate outliers in
the data.

TABLE V: Modified functions for the LoS probability distribution and performance comparisons. | Modified expression | $a$ | $b$ | MSE
---|---|---|---|---
Sigmoid | $\left(1+{a\exp(-b(\theta-a))}\right)^{-1}$ | 1.93 | 0.07 | 0.0017
Tanh | $a\frac{\exp(2b\theta)-1}{\exp(2b\theta)+1}$ | 0.9463 | 2.7972 | 0.0460
ReLu | $\max\\{0,a\theta+b\\}$ | 0.43 | 0.61 | 980.66

In addition to the modified sigmoid function suggested by authors in [22], we
propose expressions of modified tanh and ReLu functions, which are also
commonly used functions in machine learning. For a fair comparison, each of
these functions contains two undetermined parameters. Table V compares the
expressions and performances of these modified functions. The approximately
linear ReLu function is not suitable for the fitting of nonlinear curves, and
the modified sigmoid function is the most suitable one for estimating the LoS
probability distribution.

Figure 8: SCPA performance under different estimation methods.

Figure 9: A sampling of the UAV positions at different time slots for SCPA.

Fig. 9 compares the performances of SCPA under different estimation methods.
In general, more accurate LoS probability estimation leads to better
performance. The performances of the four methods are similar and
significantly better than that of the BIA.

Fig. 9 shows the positions of the UAV in different time slots obtained by SCPA
in a dynamic scenario. Assume that each user moves less than 3 meters in a
random direction (or per second), and users are always located in the outdoor
region. The results in Fig. 9 show that regardless of the used estimation
method, the positions of the UAV obtained by the SCPA do not change
significantly even after 30 slots. It suggests that only a small neighborhood
needs to be traversed by brute force.

###  V-C Performance of Algorithms Under Real-Time Search

Figure 10: User classification from coverage probability perspective when
$h=20$ m.

Figure 11: User classification from map topology perspective when $h=20$ m.

As described before, the first step of algorithms under real-time search is
user classification. Fig. 11 and Fig. 11 show the results of classification
from the coverage probability perspective and map topology perspective,
respectively. The height of the UAV is 20 m, while the degree of
classification is set to $\varepsilon=0.1$. The details of terrain-based and
non-terrain-based classification are described in subsection III-B. For the
same value of $\varepsilon$, more users will be classified into $C_{1}$ and
$C_{3}$ in terrain-based classification.

Fig. 13 provides an example of Algorithm 1 given in subsection III-C. By
executing terrain-based classification with $\varepsilon=0.1$, the two
farthest users in $C_{2}$ are selected, and the heat map of smaller average
SNR between them is given in Fig. 13. The heat map corresponds to the middle
vertical plane (the $z=0$ plane) of the two users. We can see in Fig. 13 that
the UAV starts from the LoS region on the upper right and descends in a
straight line in the direction of the steepest gradient. After arriving at the
buildings’ boundary, the UAV searches along the boundary to find the
$\gamma$-suboptimal position. Then, it enters the NLoS area of one of the
users and searches the left half branch (blue line) first and then the right
half branch (purple line) along the gradient. After confirming that the
potential LoS region is not missed, the $\gamma$-suboptimal is finally
determined as the output position.

Figure 12: An example of Algorithm 1 in a heat map.

Figure 13: Coverage probability performance comparison of all algorithms.

As shown in Fig. 13, the BIA serves as a comparison since it does not utilize
terrain information and only optimizes the horizontal coordinates. The SCPA
can further optimize the UAV deployment height based on the former algorithm
by extra terrain information collected before networking. Compared with the
above two algorithms, the MRSA further improves the coverage performance by
avoiding the shadow of buildings. The HDA has a similar coverage performance
to the MRSA, and both are close to the upper bound of coverage performance.
The upper bound is obtained by an idealized method known as brute force, which
traverses every position and leaves the position with the best overall
coverage performance.  By comparing these five algorithms, we can draw the
following conclusions:

  * •

Brute Force versus SCPA: Compared to having partial terrain information
(knowing terrain feature parameters), having complete terrain information
(precise locations of all buildings) can increase coverage probability by
about 6%.

  * •

SCPA versus BIA: Having no terrain information leads to an average loss of 9%
in coverage probability gain compared to SCPA that obtained partial terrain
information.

  * •

SCPA versus MRSA: In the case of unrestricted search length, real-time search
algorithms can gather more information during the search process than what is
contained in the feature parameters, and it brings about a 3% improvement in
coverage probability.

  * •

MRSA versus HDA: When search length is not limited, whether mastering the
values of feature parameters has little impact on real-time search.

In addition, Fig. 15 shows the boundary of the longest 20%, the mean and the
boundary of the shortest 20% of the samples of search length for the above two
algorithms. Compared to the MRSA, the HDA shortens the search distance with
the help of the collected terrain information.

Figure 14: Search lengths comparison of MRSA and HDA.

Figure 15: Influence of the number of UAVs on coverage probability performance
under multiple blockage model.

###  V-D Extended Scenarios

In this subsection, we investigate a more sophisticated and realistic scenario
in which numerous UAVs serve a region simultaneously, and the effect of
multiple blockages of the buildings is considered. The following describes the
impact of these two factors on algorithm operation.

  * •

Multiple UAVs Deployment: Firstly, the entire area is divided into sub-regions
with equal area. The number of sub-regions is equal to the number of UAVs, and
the sub-regions correspond to UAVs one by one. When the algorithm is excused
to determine the networking position, the UAV only considers the users in its
corresponding sub-region. Notice that the user will select the UAV providing
the maximum received power to be associated with, thus, the user might connect
to a UAV in another sub-region.

  * •

Modeling of Multiple Blockages: When the link between the user and the UAV is
blocked twice or more, we assume the signal from this UAV cannot be received
by the user.

Fig. 15 shows the performance of the SCPA and the MRSA while the number of
UAVs is set to one, two, and four. Compared with the results in Fig. 13, the
performance of a single UAV is slightly worse than before because of multiple
blockages, but being blocked more than once is relatively low since buildings
are sparse in our scenario. Furthermore, no matter how many UAVs are deployed,
the SCPA (blue lines) consistently outperforms the MRSA (green lines). In
addition, as the number of UAVs increases, the performance improves. This is
due to the fact that the average distance between users and UAVs decreases and
the LoS probability increases.

###  V-E Summary and Algorithm Comparison

We investigated the impact of the UAV height and the selection of mass density
function on the coverage performance of BIA. From the obtained results, we
concluded that it is more effective to reduce the weight assigned to the users
who are farther from the UAV than those who are closer. We showed the
estimation parameters $a$ and $b$ in the LoS probability model under different
regularization parameters. We recall that $a$ and $b$ are obtained by ridge
regression, which contains regularization parameters. The results have shown
that the empirical parameters do not fit well with the LoS probability
distribution.

Therefore, we tested the performance of SCPA with different values of $a$ and
$b$ under different regularization parameters, and showed that in general,
more accurate LoS probability estimation provided better, albeit minor,
improvement in the performance. We also showed that SCPA performs
significantly better than BIA. Then, we compared the coverage performance of
four algorithms with the upper bound, that is, the brute force method. The BIA
serves as a lower bound since it does not utilize terrain information and only
optimizes the horizontal coordinates. On the other hand, the SCPA can further
optimize the UAV deployment height based on the former algorithm by extra
terrain information collected before networking, and therefore, has better
performance than the BIA. Compared with the above two algorithms, the MRSA
further improves the coverage performance by avoiding the shadow of buildings.
The HDA has a similar coverage performance to the MRSA, and both are close to
the upper bound of coverage performance. We compared the search lengths of HDA
and MRSA since they are two real-time search methods, the difference between
them is that the former has offline terrain information and the latter does
not. Compared to the MRSA, the HDA shortens the search distance with the help
of the collected terrain information.

Finally, we extended the scenario using only one UAV to a general scenario
using multiple UAVs. We also extended the original blockage model from using
only LoS and NLoS into a multiple blockage model, which includes no blockage,
blockage once, and blockage twice or more than twice. We showed that as the
number of UAVs increases, the performance improves. However, the performance
of a single UAV is slightly worse compared to the scenario using only LoS and
NLoS. This minor decrease in performance is caused by the multiple blockages,
however, the probability of being blocked more than once is relatively low
since buildings are sparse in the considered scenario.

##  VI Conclusion

In this paper, two methods were proposed to improve the coverage of the UAV
network in combination with the actual terrain information before and during
networking. Before networking, the UAV can estimate the probability that UAV
and users are blocked by buildings at different elevation angles through
offline information collection. Furthermore, we derived the analytical
expression of coverage probability and classified users accordingly.
Classification results can simplify the optimization problem and narrow the
range of users that UAVs need to focus on. During networking, we analyzed the
optimality of UAV networking position in simple scenes, and further designs
the real-time search method. Furthermore, we designed four algorithms
according to whether the above two UAV deployment methods are used or not.
Finally, in order to obtain more feasible experimental results, multiple UAV
deployments and multiple blockages are considered as further extensions.

##  Appendix A Proof of Theorem 1

According to the definition in (3), the coverage probability for a user whose
distance from the UAV is $r_{k}$ can be expressed as,

$\begin{split}&P_{C,k}\left(h,r_{k}\right)=P_{\rm{LoS}}\left(h,r_{k}\right)\mathbb{P}\left[\frac{S_{{\rm{LoS}},k}\left(r_{k}\right)}{\sigma^{2}}>\gamma\right]\\\
&+P_{\rm{NLoS}}\left(h,r_{k}\right)\mathbb{P}\left[\frac{S_{{\rm{NLoS}},k}\left(r_{k}\right)}{\sigma^{2}}>\gamma\right]\\\
&\overset{(a)}{=}P_{\rm{LoS}}\left(h,r_{k}\right)\frac{\Gamma_{u}\left(m_{\rm{LoS}},m_{\rm{LoS}}\eta_{\rm{LoS}}^{-1}\zeta^{-1}\gamma\,r_{k}^{\alpha_{\rm{LoS}}}\sigma^{2}\right)}{\Gamma\left(m_{\rm{LoS}}\right)}\\\
&+P_{\rm{NLoS}}\left(h,r_{k}\right)\frac{\Gamma_{u}\left(m_{\rm{NLoS}},m_{\rm{NLoS}}\eta_{\rm{NLoS}}^{-1}\zeta^{-1}\gamma\,r_{k}^{\alpha_{\rm{NLoS}}}\sigma^{2}\right)}{\Gamma\left(m_{\rm{NLoS}}\right)}\\\
&\overset{(b)}{=}P_{\rm{LoS}}\left(h,r_{k}\right)\exp\left(-\mu_{\rm{LoS}}\left(r_{k}\right)\right)\sum_{n=0}^{m_{\rm{LoS}}-1}\frac{\mu_{\rm{LoS}}\left(r_{k}\right)^{n}}{n!}\\\
&+P_{\rm{NLoS}}\left(h,r_{k}\right)\exp\left(-\mu_{\rm{NLoS}}\left(r_{k}\right)\right)\sum_{n=0}^{m_{\rm{NLoS}}-1}\frac{\mu_{\rm{NLoS}}\left(r_{k}\right)^{n}}{n!},\end{split}$
(12)

where $(a)$ is from the complementary cumulative distribution function (CCDF)
of the Gamma function
$\overline{F}_{G}\left(g\right)=\frac{\Gamma_{u}\left(m,mg\right)}{\Gamma\left(m\right)}$,
and $\Gamma_{u}\left(m,mg\right)=\int_{mg}^{\infty}t^{m-1}e^{-t}dt$ is the
upper incomplete Gamma function. $(b)$ is derived by the definition
$\frac{\Gamma_{u}\left(m,mg\right)}{\Gamma\left(m\right)}=\exp(-g)\sum_{n=0}^{m-1}\frac{g^{n}}{n!}$
and substituting
$\mu_{Q}\left(r_{k}\right)=\eta_{Q}^{-1}\zeta^{-1}\gamma\,r_{k}^{\alpha_{Q}}\sigma^{2},\
Q=\left\\{\rm{LoS},\rm{NLoS}\right\\}$ into the expression in $(a)$.

##  Appendix B Proof of Proposition 1

In scenario I, the height of the optimal position is $h_{\min}$ since there
are no buildings. Therefore, only the horizontal position of the UAV remains
to be optimized.

When a UAV serves only one user, the optimal horizontal position is the user’s
horizontal position. In other words, the UAV hovers over the user.

When covering two users who are close to each other, the optimal horizontal
position of the UAV might be the midpoint of the two users. When two users are
far apart, there are two optimal locations where there is symmetry about the
midpoint of the two users. When two users are infinitely far apart, the UAV
hovers over one of them. In either case, the optimal location is on the line
segment connecting the two users and is therefore linearly optimal.

When the UAV provides coverage to more than three users, we simplify the
optimization objectives as follows:

$\begin{split}&P_{C}\left(h_{\min};r_{1},...,r_{K}\right)=\frac{1}{K}\sum_{k=1}^{K}\mathbb{P}\left[\frac{S_{Q,k}\left(r_{k}\right)}{\sigma^{2}}>\gamma\right]\\\
&\overset{(a)}{\approx}\frac{1}{K}\sum_{k=1}^{K}v\sqrt{r_{k}^{2}+h_{\min}^{2}}+w\overset{(b)}{\approx}\frac{1}{K}\sum_{k=1}^{K}vr_{k}+w,\end{split}$
(13)

where $(a)$ is obtained by approximating the coverage probability as a linear
function of the Euclidean distance between the UAV and the users. Through
classification, the coverage probabilities of users in $C_{2}$ locate in the
approximately linear region of the covering probability function, so this
assumption is reasonable. $(b)$ is based on the fact that $r_{k}$ is often
much larger than $h_{\min}$ for $C_{2}$ users.

When $K=3$, the optimal position is called the Fermat-Weber point. When $K\geq
4$, the authors in [40] proposed a method to approach the optimal position
based on the above approximation. Since the number of iterations is finite,
the search trajectory is linear. However, it requires that the polygon
connected by users is a convex polygon. In addition, this linear search
trajectory cannot approach with an arbitrarily small precision, therefore, it
is not linear optimal.

##  Appendix C Proof of Proposition 2

The formula (1) indicates that for a given UAV position, the set of positions
on the ground satisfying $\overline{\rm{SNR}}_{Q}(r)\geq\gamma$ is a circle.
Therefore, if all the users are located inside a circle or on its boundary,
then the center of that circle is one of the suboptimal solutions. This
problem is called the circle covering problem. In [42], the author provides a
minimum circle covering algorithm. The algorithm keeps two points (users) on
the boundary of the covering circle, repeatedly shrink the circle and replaces
the boundary points (users) until the optimal covering circle is found.
Through this algorithm, the UAV can find the suboptimal location without
searching, thus, a $\gamma$-suboptimal strategy $\mathscr{L}$ exists.

##  Appendix D Proof of Proposition 3

Before starting the proof, we denote the UAV location that is not blocked by
the building for both users as LoS position, otherwise NLoS position. Under
the condition that $2\rho^{*}<d$, the proof of Proposition 3 is divided into
the following four steps:

  * •

Step 1: The $\gamma$-suboptimal position is located on the boundaries caused
by the shadows of the buildings, or on the plane of height $h_{\min}$.

  * •

Step 2: In most cases, the $\gamma$-suboptimal position is located in the
plane $z=0$.

  * •

Step 3: The output $(\rho^{*},\theta^{*},0)$ is the $\gamma$-suboptimal
position.

  * •

Step 4: The search length of Algorithm 1 is linear with respect to the height
of UAV $h$ and the distance between users $d$.

Step 1: For LoS positions, the average SNR of both users is enhanced with the
decrease of the UAV height. For NLoS positions, the average SNR of the users
with weak signal power is enhanced with the decrease of the UAV height. Due to
the minimum height restrictions, the $\gamma$-suboptimal position is located
on the boundaries caused by the shadows of the buildings (corresponding to the
LoS positions), or on the plane of height $h_{\min}$ (corresponding to the
NLoS positions).

Step 2: The right side of the Fig. 2 shows an inclined plane that can be
uniquely determined by three points on the plane, that is, the coordinates of
two users $\left(0,0,\pm\frac{d}{2}\right)$ and an arbitrary point
$\left(\rho_{0},\theta_{0},0\right)$, where $\theta_{0}$ is the inclination
between the plane and the ground. Using basic geometry, we can draw the
following conclusion: for an arbitrarily small $\Delta\rho\rightarrow 0$, if
$\left(\rho_{0}^{*}+\Delta\rho,\theta_{0},0\right)$ is a LoS position and
$\left(\rho_{0}^{*}-\Delta\rho,\theta_{0},0\right)$ is a NLoS position
($\rho_{0}^{*}\sin\theta_{0}>h_{\min}$), the smaller average SNR of any LoS
position in the plane with inclination of $\theta_{0}$ is less than
$\overline{\rm{SNR}}_{\rm{LoS}}\left(\max\left\\{\frac{2\rho_{0}^{*}d}{\sqrt{\left(2\rho_{0}^{*}\right)^{2}+d^{2}}},\sqrt{\left(\rho_{0}^{*}\right)^{2}+\frac{d^{2}}{4}}\right\\}\right)$.
When $2\rho_{0}^{*}<d$, $\left(\rho_{0}^{*},\theta_{0},0\right)$ is the
position with the maximum smaller average SNR, denoted as the inclined-
suboptimal position. Because in most cases, $2\rho_{0}^{*}<d$ is satisfied (in
simulation, $d\approx 120$m and $\rho^{*}\approx 20$m usually), thus, the
suboptimal position is located in the plane $z=0$.

Step 3: We extend the trajectory to the left and right of the straight line
$h=h_{\min}$ in the plane $z=0$. According to the monotonically decreasing
value of $\rho$ in algorithm 1, $(\rho^{*},\theta^{*},0)$ has the shortest
distance to the users among the points on the search trajectory. It is
straightforward to prove that any point on the extended trajectory of
Algorithm 1 is less than or equal to the distance from the corresponding
inclined-suboptimal position to the user. Since $(\rho^{*},\theta^{*},0)$ is a
LoS position, and has the shortest distance to the users among the inclined-
suboptimal positions, $(\rho^{*},\theta^{*},0)$ is the $\gamma$-suboptimal
position.

Step 4: Since one lower bound of the search trajectory is
$2\arccos\left(\frac{h_{\min}}{\rho^{*}}\right)\sqrt{\left(\rho^{*}\right)^{2}+\frac{d^{2}}{4}}$
and one upper bound is
$2\arccos\left(\frac{h_{\min}}{\rho_{\max}}\right)\sqrt{\left(\rho_{\max}\right)^{2}+\frac{d^{2}}{4}}$,
where $\rho_{\max}$ is $\rho_{u}$ after step (6) in Algorithm 1, thus the
trajectory is linear. Finally, note that although we have shown that the
$\gamma$-subpotimal position might be located on the boundary and in the plane
$z=0$ in most cases, traversing along the boundary does not guarantee
linearity. Because boundaries can be complex and finding an upper bound of the
search trajectory is not guaranteed.

## References

  * [1] Z. Hua, Y. Lu, G. Pan, K. Gao, D. B. d. Costa, and S. Chen, “Computer vision-aided mmWave UAV communication systems,” _IEEE Internet Things J._ , vol. 10, no. 14, pp. 12 548–12 561, Mar. 2023.
  * [2] B. Shang, E. S. Bentley, and L. Liu, “UAV swarm-enabled aerial reconfigurable intelligent surface: Modeling, analysis, and optimization,” _IEEE Trans. Commun._ , vol. 71, no. 6, pp. 3621–3636, May 2022.
  * [3] Z. Lou, B. E. Y. Belmekki, and M.-S. Alouini, “HAPS in the non-terrestrial network nexus: Prospective architectures and performance insights,” _arXiv preprint arXiv:2310.09659_ , 2023.
  * [4] R. Wang, M. A. Kishk, and M.-S. Alouini, “Ultra-dense LEO satellite-based communication systems: A novel modeling technique,” _IEEE Commun. Mag._ , vol. 60, no. 4, pp. 25–31, Apr. 2022.
  * [5] Q. Huang, B. E. Y. Belmekki, A. M. Eltawil, and M.-S. Alouini, “System-level metrics for non-terrestrial networks under stochastic geometry framework,” _arXiv preprint arXiv:2302.03376_ , 2023.
  * [6] R. Wang, M. A. Kishk, and M.-S. Alouini, “Reliability analysis of multi-hop routing in multi-tier leo satellite networks,” _arXiv preprint arXiv:2303.02286_ , 2023.
  * [7] M. A. Imran, A. Turkmen, M. Ozturk, J. Nadas, and Q. H. Abbasi, “Seamless indoor/outdoor coverage in 5G,” _Wiley 5G Ref: The Essential 5G Reference Online_ , pp. 1–23, May 2020.
  * [8] J. Du, D. Chizhik, R. A. Valenzuela, R. Feick, G. Castro, M. Rodriguez, T. Chen, M. Kohli, and G. Zussman, “Directional measurements in urban street canyons from macro rooftop sites at 28 GHz for 90% outdoor coverage,” _IEEE Trans. Antennas Propag._ , vol. 69, no. 6, pp. 3459–3469, Dec 2020.
  * [9] R. Wang, M. A. Kishk, and M.-S. Alouini, “Stochastic geometry-based low latency routing in massive LEO satellite networks,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 58, no. 5, pp. 3881–3894, Aug. 2022.
  * [10] K. B. Letaief, W. Chen, Y. Shi, J. Zhang, and Y.-J. A. Zhang, “The roadmap to 6G: AI empowered wireless networks,” _IEEE Commun. Mag._ , vol. 57, no. 8, pp. 84–90, Aug. 2019.
  * [11] B. E. Y. Belmekki and M.-S. Alouini, “Unleashing the potential of networked tethered flying platforms: Prospects, challenges, and applications,” _IEEE Open J. Veh. Technol._ , vol. 3, pp. 278–320, May 2022.
  * [12] R. Wang, M. A. Kishk, and M.-S. Alouini, “Resident population density-inspired deployment of K-tier aerial cellular network,” _IEEE Trans. Wirel. Commun._ , vol. 22, no. 11, pp. 7989–8002, Mar. 2023.
  * [13] E. Yaacoub and M.-S. Alouini, “A key 6G challenge and opportunity-connecting the base of the pyramid: A survey on rural connectivity,” _Proceedings of the IEEE_ , vol. 108, no. 4, pp. 533–582, Mar. 2020.
  * [14] R. Wang, M. A. Kishk, and M.-S. Alouini, “Evaluating the accuracy of stochastic geometry based models for LEO satellite networks analysis,” _IEEE Commun. Lett._ , vol. 26, no. 10, pp. 2440–2444, Jul. 2022.
  * [15] A. Elzanaty, A. Guerra, F. Guidi, D. Dardari, and M.-S. Alouini, “Toward 6G holographic localization: Enabling technologies and perspectives,” _IEEE Internet Things Mag._ , vol. 6, no. 3, pp. 138–143, Apr. 2023.
  * [16] Z. Lou, B. E. Y. Belmekki, and M.-S. Alouini, “Coverage analysis of hybrid RF/THz networks with best relay selection,” _IEEE Commun. Lett._ , vol. 27, no. 6, pp. 1634–1638, Apr. 2023.
  * [17] A. A. Khuwaja, Y. Chen, N. Zhao, M.-S. Alouini, and P. Dobbins, “A survey of channel modeling for UAV communications,” _IEEE Commun. Surv. Tutor._ , vol. 20, no. 4, pp. 2804–2821, Jul. 2018.
  * [18] A. A. Zaid, B. E. Y. Belmekki, and M.-S. Alouini, “eVTOL communications and networking in UAM: Requirements, key enablers, and challenges,” _IEEE Commun. Mag._ , vol. 61, no. 8, pp. 154–160, Aug. 2023.
  * [19] Z. Lou, A. Elzanaty, and M.-S. Alouini, “Green tethered UAVs for EMF-aware cellular networks,” _IEEE Trans. Green Commun. Netw._ , vol. 5, no. 4, pp. 1697–1711, Aug. 2021.
  * [20] A. A. Zaid, B. E. Y. Belmekki, and M.-S. Alouini, “Aerial-aided mmWave VANETs using NOMA: Performance analysis, comparison, and insights,” _IEEE Trans. Veh. Technol._ , pp. 1–16, Nov. 2023.
  * [21] R. Wang, W. U. Mondal, M. A. Kishk, V. Aggarwal, and M.-S. Alouini, “Terrain-based coverage manifold estimation: Machine learning, stochastic geometry, or simulation?” _arXiv preprint arXiv:2312.01826_ , 2023.
  * [22] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” _IEEE Wirel. Commun. Lett._ , vol. 3, no. 6, pp. 569–572, Jul. 2014.
  * [23] D. Gesbert, O. Esrafilian, J. Chen, R. Gangula, and U. Mitra, “UAV-aided RF mapping for sensing and connectivity in wireless networks,” _IEEE Wirel. Commun._ , vol. 30, no. 4, pp. 116–122, Aug. 2023.
  * [24] J. Chen and D. Gesbert, “Optimal positioning of flying relays for wireless networks: A LOS map approach,” in _Proc. IEEE Int. Conf. Commun. (ICC), Paris, France_ , Jul. 2017, pp. 1–6.
  * [25] H. Nawaz, H. M. Ali, and A. A. Laghari, “UAV communication networks issues: A review,” _Arch. Comput. Methods Eng._ , vol. 28, no. 3, pp. 1349–1369, Mar. 2021.
  * [26] S. Dang, O. Amin, B. Shihada, and M.-S. Alouini, “What should 6G be?” _Nat. Electron._ , vol. 3, no. 1, pp. 20–29, Jan. 2020.
  * [27] Q. Feng, E. K. Tameh, A. R. Nix, and J. McGeehan, “WLCp2-06: Modelling the likelihood of line-of-sight for air-to-ground radio propagation in urban environments,” in _Proc. IEEE Global Comm. Conf. (Globecom), San Francisco, CA, USA_ , Apr. 2006, pp. 1–5.
  * [28] A. Al-Hourani, S. Kandeepan, and A. Jamalipour, “Modeling air-to-ground path loss for low altitude platforms in urban environments,” in _Proc. IEEE Global Comm. Conf. (Globecom), Austin, TX, USA_ , Dec. 2014, pp. 2898–2904.
  * [29] M. Matracia, M. A. Kishk, and M.-S. Alouini, “Coverage analysis for UAV-assisted cellular networks in rural areas,” _IEEE Open J. Veh. Technol._ , vol. 2, pp. 194–206, Apr. 2021.
  * [30] S. Yuan, Y. Sun, and M. Peng, “Joint network function placement and routing optimization in dynamic software-defined satellite-terrestrial integrated networks,” _IEEE Trans. Wirel. Commun._ , pp. 1–1, Oct. 2023.
  * [31] 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 Wirel. Commun. Lett._ , vol. 6, no. 4, pp. 434–437, May 2017.
  * [32] X. He, W. Yu, H. Xu, J. Lin, X. Yang, C. Lu, and X. Fu, “Towards 3D deployment of UAV base stations in uneven terrain,” in _Proc. IEEE 27th International Conference on Computer Communication and Networks (ICCCN), Hangzhou, China_ , Oct. 2018, pp. 1–9.
  * [33] J. Chen and D. Gesbert, “Efficient local map search algorithms for the placement of flying relays,” _IEEE Trans. Wirel. Commun._ , vol. 19, no. 2, pp. 1305–1319, Nov. 2020.
  * [34] J. Chen, U. Mitra, and D. Gesbert, “3D urban UAV relay placement: Linear complexity algorithm and analysis,” _IEEE Trans. Wirel. Commun._ , vol. 20, no. 8, pp. 5243–5257, Mar. 2021.
  * [35] ——, “3D urban UAV relay placement: Linear complexity algorithm and analysis,” _IEEE Trans. Wirel. Commun._ , vol. 20, no. 8, pp. 5243–5257, Mar. 2021.
  * [36] M. Alzenad and H. Yanikomeroglu, “Coverage and rate analysis for vertical heterogeneous networks (VHetNets),” _IEEE Trans. Wirel. Commun._ , vol. 18, no. 12, pp. 5643–5657, Sep. 2019.
  * [37] D. Wackerly, W. Mendenhall, and R. L. Scheaffer, _Mathematical Statistics with Applications_. Cengage Learning, 2014.
  * [38] R. I. Bor-Yaliniz, A. El-Keyi, and H. Yanikomeroglu, “Efficient 3-D placement of an aerial base station in next generation cellular networks,” in _Proc. IEEE Int. Conf. Commun. (ICC), Kuala Lumpur, Malaysia_ , May 2016, pp. 1–5.
  * [39] Y. Zou and K. Chakrabarty, “A distributed coverage-and connectivity-centric technique for selecting active nodes in wireless sensor networks,” _IEEE Trans. Comput._ , vol. 54, no. 8, pp. 978–991, Jun. 2005.
  * [40] P. Carmi, S. Har-Peled, and M. J. Katz, “On the Fermat–Weber center of a convex object,” _Computational Geometry_ , vol. 32, no. 3, pp. 188–195, 2005.
  * [41] T. Zhang, J. Zhao, L. An, and D. Liu, “Energy efficiency of base station deployment in ultra dense HetNets: A stochastic geometry analysis,” _IEEE Wirel. Commun. Lett._ , vol. 5, no. 2, pp. 184–187, Jan. 2016.
  * [42] D. W. Hearn and J. Vijay, “Efficient algorithms for the (weighted) minimum circle problem,” _Operations Research_ , vol. 30, no. 4, pp. 777–795, Aug. 1982.