# Energy-Efficient URLLC Service Provision via a Near-Space Information
Network

Puguang An, Xianbin Cao, , Peng Yang, , Kun Guo, , Yue Gao, , and Tony Q. S.
Quek,  P. An is with the School of Electronic and Information Engineering,
Beihang University, Beijing 100191, China. X. Cao and P. Yang are with the
School of Electronic and Information Engineering, Beihang University, Beijing
100191, China, also with the Key Laboratory of Advanced Technology of Near
Space Information System (Beihang University), Ministry of Industry and
Information Technology of China, Beijing 100191, China, and also with the Peng
Cheng Laboratory, Department of Mathematics and Theories, Shenzhen, Guangdong
518055, China. K. Guo is with the School of Communication and Electronic
Engineering, East China Normal University, Shanghai 200062, China. Y. Gao is
with the School of Computer Science, Fudan University, Shanghai 200433, China.
T. Q. S. Quek is with the Information Systems Technology and Design, Singapore
University of Technology and Design, 487372 Singapore. This work has been
submitted to the IEEE for possible publication. Copyright may be transferred
without notice, after which this version may no longer be accessible.

###### Abstract

The integration of a near-space information network (NSIN) with the
reconfigurable intelligent surface (RIS) is envisioned to significantly
enhance the communication performance of future wireless communication systems
by proactively altering wireless channels. This paper investigates the problem
of deploying a RIS-integrated NSIN to provide energy-efficient, ultra-reliable
and low-latency communications (URLLC) services. We mathematically formulate
this problem as a resource optimization problem, aiming to maximize the
effective throughput and minimize the system power consumption, subject to
URLLC and physical resource constraints. The formulated problem is challenging
in terms of accurate channel estimation, RIS phase alignment, theoretical
analysis, and effective solution. We propose a joint resource allocation
algorithm to handle these challenges. In this algorithm, we develop an
accurate channel estimation approach by exploring message passing and optimize
phase shifts of RIS reflecting elements to further increase the channel gain.
Besides, we derive an analysis-friend expression of decoding error probability
and decompose the problem into two-layered optimization problems by analyzing
the monotonicity, which makes the formulated problem analytically tractable.
Extensive simulations have been conducted to verify the performance of the
proposed algorithm. Simulation results show that the proposed algorithm can
achieve outstanding channel estimation performance and is more energy-
efficient than diverse benchmark algorithms.

###### Index Terms:

Near-space information network, reconfigurable intelligent surface,
electromagnetic channel estimation, URLLC, energy efficiency

##  I Introduction

The upcoming sixth generation (6G) communication systems have the
unprecedented requirements for ultra-reliability, ultra-low latency, and
extremely high data rates that the fifth generation (5G) communication systems
cannot well fulfill [1]. Ultra-reliable and low-latency communications
(URLLC), which is one of the three pillar services for the 5G systems, will
still act as a pillar service for the 6G systems. URLLC can empower diverse
mission-critical applications including intelligent transportation systems,
autonomous driving, telemedicine, industrial automation, Tactile Internet,
real-time control in Digital Twins, and Metaverse.

It’s widely considered that space-air-ground integrated communication systems
are a crucial system architecture that can satisfy the unprecedented
requirements of 6G systems. A near-space information network (NSIN), composed
of various near-space platforms (e.g., airships), medium/low-altitude
platforms (e.g., medium/low-altitude unmanned aerial vehicles (UAVs)), is a
network system that can acquire, transmit, and process space electromagnetic
(EM) signals in real-time. It’s the core component of the space-air-ground
integrated communication systems. Compared to space-based networks, NSIN is
much closer to terrestrial users and can provide high-capacity communication
services. Compared with terrestrial networks, NSIN can quickly extend the
communication distance and is robust to unexpected events (e.g., diverse
natural disasters). Thus, NSIN can provide wide-area and reliable
communication coverage. In this regard, NSIN integrates the advantages of
space-based networks and terrestrial networks and becomes an important
supplement to terrestrial networks and space-based networks. It’s envisioned
that resorting to NSIN will become a new regime for enabling URLLC [2, 3].

Actually, the authors have demonstrated the feasibility of delivering reliable
and low-latency control-and-non-payload EM signals using one of our
manufactured near-space airships during the execution of many projects on
airship-supported intelligent transportation systems and meteorological
observations over the past few years.

Reconfigurable intelligent surface (RIS) has emerged as a promising paradigm
to achieve smart wireless reflection channels for future wireless
communication systems. RIS is a planar surface consisting of many passive
reflecting elements for EM waves, each of which can separately change the
amplitude and/or phase of impinging EM signals. With the deployment of RIS,
wireless channels between transmitters and receivers can be flexibly
reconfigured to achieve desired realizations and distributions. Thus, RIS
provides a new paradigm to fundamentally address the complex EM signal
interference and channel fading issues and potentially achieves a significant
improvement in terms of throughput and reliability of wireless communication
systems [4].

Different from existing advanced communication techniques (e.g., beamforming,
dynamic rate control) to realize ultra-reliable and high-capacity transmission
by adapting to channel conditions, both NSIN and RIS can proactively alter the
wireless channels to significantly improve communication performance.
Specifically, NSIN can achieve rich LoS transmission through
trajectory/location control. RIS can create desired LoS links for EM wave
propagation and eliminate EM signal interference via smart EM signal
reflection.

The integration of NSIN with RIS achieves complementary advantages, which will
significantly increase the design flexibility of NSIN, reduce cost, and
enhance communication performance. Through mounting RIS on an airborne
platform of NSIN, RIS can achieve a 360° panoramic full-angle reflection
toward target nodes. By deploying RIS-integrated NSIN, the shortcoming of
wide-area interference inherent in NSIN can be effectively compensated for,
and the spectral efficiency of NSIN can be improved. Meanwhile, the
integration of NSIN with RIS will cost-effectively improve the energy
efficiency of NSIN due to the low hardware and energy costs of RIS.

###  I-A Related work

Presently, the research on the URLLC-enabled UAV network or NSIN has attracted
much attention, as outlined in [3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. For
example, the work in [3] proposed an NSIN with low/medium/high altitude
platforms to support URLLC applications. The author explored a complex network
optimization tool, i.e., polychromatic sets, to analyze the performance of
each link. Based on the analysis results, they designed distributed and
centralized link selection schemes to select end-to-end links that could
guarantee the stringent latency and reliability requirements of URLLC. The
work in [5] considered a UAV relay to communicate short URLLC control packets
between a controller and many mobile robots. In particular, the author studied
the optimization of UAV deployment location, beamwidth, transmit power, and
blocklength allocated to each robot and proposed an iterative optimization
method to obtain the optimal solution. The work in [6] investigated the
optimization problem of UAV deployment altitude, blocklength, transmission
power, and update generation rate to guarantee the freshness of received
information at the ground destination node in an URLLC-enabled UAV network.
The work in [7] analyzed channel quality indicator (CQI) aging effects in
URLLC under the finite blocklength regime and utilized a Recurrent Neural
Network (RNN) approach to predict the next CQI for UAV control information
exchange. The author in [8] studied the average packet error probability and
effective throughput of URLLC control links in UAV communications. They
derived the closed-form expression of packet error probability and throughput
using Gaussian-Chebyshev quadrature and obtained the optimal packet length
using an exhaustive search. Besides, observing the high possibility of
establishing LoS links between UAVs and ground users, the author in [9]
investigated the potential of using UAVs in supporting URLLC. They
characterized the transmission latency, reliability, and network availability
of UAV networks, optimized the UAV deployment altitude, and allocated
bandwidth that minimized the required total bandwidth for URLLC.

However, the work in [3, 5, 6, 7, 8, 9, 10, 11, 12] did not investigate the
NSIN- or UAV-assisted URLLC systems with the aid of RIS. RIS is an important
emerging technique that can provide optimized energy and spectral efficiency
for UAV-assisted URLLC systems in a cost-effective manner. Recently, many
researchers have proposed to design URLLC systems by the joint optimization of
passive beamforming and UAV resources [15, 16, 17, 18]. For instance, the
authors in [15] proposed a UAV-RIS system to support the stringent URLLC
constraints and formulated an optimization framework in terms of UAVs’
locations, the phase shift of RIS, and the blocklength of URLLC. To solve the
formulated problem, a deep neural network (DNN) architecture was designed to
obtain 3-D deployment locations of UAVs, and an iterative optimization scheme
was explored to optimize the phase shift and blocklength. The work in [16]
considered a scenario of deploying a UAV and a ground RIS to deliver URLLC
packets between ground Internet-of-Things (IoT) devices. The authors proposed
to jointly optimize the passive beamforming, the deployment location of a UAV,
and the blocklength to minimize the packet decoding error probability (DEP). A
computationally efficient Nelder–Mead simplex approach was utilized to solve
the joint optimization problem. The authors in [17] developed a joint passive
beamforming and UAV trajectory optimization approach to support URLLC
applications over UAV/RIS-based systems. They first built UAV/RIS-based
channel models in a finite blocklength regime. Then, they formulated a finite-
blocklength coding-based energy-efficiency maximization problem. Next, an
alternative optimization scheme was utilized to optimize UAV transmit power,
passive beamforming, and UAV trajectory alternatively. Additionally, the
authors in [18] investigated the benefits of UAV/RIS-based systems to
passively relay low-latency information sampled by Internet of Things (IoT)
devices to a base station (BS). For this aim, they formulated a joint
optimization problem of UAV altitude, the phase shift of RIS, and
communication scheduling, aiming to minimize the expected sum Age-of-
Information (AoI) under a reliable decoding constraint. They also designed a
DRL framework based on proximal policy optimization to solve the optimization
problem.

###  I-B Motivations and Contributions

Although the issue of designing UAV/RIS-based systems to support URLLC
applications was investigated in the works [15, 16, 17, 18], they did not
consider the scenario of integrating high altitude platforms (HAPs), UAVs, and
aerial RIS to provide wide-area and robust communication coverage in an
energy- and cost-efficient manner. Most recently, the works in [19, 20]
investigated the integration of NSIN and RIS. For example, the authors in [19]
considered a novel wireless architecture composed of a RIS-enabled HAP and
multiple UAVs. The HAP-RIS could overcome the limitations of conventional
terrestrial RIS and achieve rich LoS and full-area coverage. The UAVs that
would receive backhaul reflection signals from HAP-RIS were deployed to serve
terrestrial flash crowd traffic. They formulated and solved a joint
optimization problem of HAP deployment location and phase shift of RIS to
increase the energy efficiency of UAVs. The work in [20] considered a scenario
of deploying a UAV-RIS to assist the HAP downlink transmission. The authors
formulated an optimization problem of maximizing the signal-to-noise ratio
(SNR) of ground users by jointly optimizing the UAV trajectory and the phase
shift of UAV-RIS. They developed a model-free reinforcement learning framework
to learn the UAV trajectory and adjust the phase shift according to the
learned UAV trajectory.

The works [19, 20], however, investigated the high throughput scenario instead
of the URLLC scenario. Moreover, they did not studied the crucial HAP-UAV
channel estimation problem. A LoS propagation assumption between a HAP and a
low-altitude UAV was enforced in [19, 20]. This assumption will be
unreasonable as the received signals by a low-altitude UAV include both LoS
and non-line-of-sight (NLoS) (e.g., reflected/scattered signals) components.
Therefore, it’s essential to study the HAP-UAV channel estimation issue.

To tackle the above issues, this paper proposes to design a RIS-integrated
NSIN to support energy-efficient URLLC services, which leverages the
advantages of power control and RIS. The main contributions of this paper are
summarized as follows:

  * •

We formulate the problem of providing energy-efficient URLLC services using a
RIS-integrated NSIN as a resource optimization problem. The goal of this
problem is to maximize the effective throughput and minimize the system power
consumption, subject to URLLC and physical resource constraints. The
formulated problem is challenging in terms of channel estimation, RIS phase
shift optimization, theoretical analysis, and effective solution.

  * •

To handle these challenges, we propose a joint resource allocation algorithm
in which we first develop a probabilistic channel model to capture the
sparsity of the HAP-UAV channel and then design a novel channel estimation
approach by message passing to accurately estimate the HAP-UAV channel.

  * •

Second, we apply a maximum-ratio transmission (MRT) strategy to maximize the
signal-to-noise ratio (SNR) and design a RIS phase shift optimization
strategy, resulting in a suboptimal HAP-UAV channel gain.

  * •

Third, we derive an analysis-friend expression of DEP by Taylor expansion to
make the formulated problem analytically tractable and help to lower the DEP
value.

  * •

Fourth, through analyzing the monotonicity of the objective function in the
feasible region, we decouple the decision variables and horizontally decompose
the problem into a BS-Layer optimization problem and a UAV-Layer optimization
problem. An iterative optimization strategy is then explored to optimize
transmit power and blocklength alternately.

  * •

Finally, the proposed algorithm is compared with diverse benchmarks to verify
its performance, and the impact of various design parameters is
comprehensively discussed. Simulation results show that the proposed algorithm
can achieve accurate HAP-UAV channel estimation and can provide energy-
efficient URLLC services.

##  II System model and problem formulation

###  II-A System Model

We consider a downlink transmission scenario based on a RIS-integrated NSIN
communication system for energy-efficient URLLC service provision. As
illustrated in Fig. 1, the system mainly includes a base station (BS), a RIS-
mounted HAP, a low-attitude UAV, and multiple ground remote URLLC robots. The
BS is designated to transmit ultra-reliable and low-latency signals towards
these remote robots for accomplishing some critical tasks (e.g., grid,
pipeline, security, and remote railway patrol, farming, and mining).

Figure 1: A downlink URLLC transmission scenario with a RIS-integrated NSIN
communication system.

A direct propagation path between the BS and a remote URLLC robot, however,
will be blocked by numerous physical obstructions on the ground (e.g.,
skyscrapers, mountains, and trees) over the remote propagation path. Then, the
RIS-mounted HAP and the UAV can work cooperatively to establish ultra-reliable
near-space communication links between the BS and multiple remote URLLC
robots. The RIS-mounted HAP acts as a passive flying relay, and the UAV acts
as an active flying relay. In this system, the BS is located at the origin and
equipped with a uniform linear-array (ULA) antenna with $M$ antenna array
elements. The antenna gain equals $G$, and the spacing of adjacent antenna
elements is $d_{\rm BS}$. There are $K$ URLLC robots with a single antenna
configuration, and they are randomly and uniformly distributed in a considered
area $\mathcal{G}$. The URLLC robot set is denoted as
${\mathcal{K}}=\left\\{1,...,K\right\\}$, and denote the coordinate of the
$k$-th robot by ${\bm{p}}_{G}=[x_{G,k},y_{G,k}]^{\rm T}$. The UAV is equipped
with an omnidirectional antenna, hovering in the air with a fixed three-
dimensional (3-D) coordinate ${\bm{p}}_{U}=[x_{U},y_{U},h_{U}]^{\rm T}$. A
decoding-and-forward (DF) relay strategy is adopted to flexibly control the
blocklength for the BS-HAP-UAV link and the UAV-to-ground (UtG) link.
Moreover, the frequency division multiple access (FDMA) manner is exploited to
eliminate intra-cell interference.

For the RIS-mounted HAP, we assume that a ULA-structured RIS is deployed
directly below the HAP to provide wide-area and full-angle coverage. The RIS
includes $N$ passive reflecting elements, separated by the spacing $d_{\rm
RIS}$. We assume that the BS array is deployed on the origin. The HAP hovers
at a fixed point ${\bm{p}}_{H}=[x_{H},y_{H},h_{H}]^{\rm T}$ in the
stratosphere. We assume that all elements of the RIS are turned on without
power amplification. Thus, the diagonal phase shift matrix of the RIS can be
modeled as [21, 22, 23]

${\rm\bm{\Theta}}={\rm
diag}\left\\{e^{j\theta_{1}},e^{j\theta_{2}},...,e^{j\theta_{N}}\right\\}\in\mathbb{C}^{N\times
N}$ (1)

where ${\theta_{n}}\in(0,2\pi]$ is the phase shift parameter of the $n$-th
element. During the first transmission phase, the BS will transmit URLLC
signals to the UAV with the assistance of the HAP relay. In this phase, a
cascading channel comprised of the BS-HAP channel and the HAP-UAV channel is
constructed. Denote the BS-HAP channel by ${\bm{H}}\in\mathbb{C}^{N\times M}$,
and the HAP-UAV channel by ${\bm{h}}\in\mathbb{C}^{N}$. Therefore, the
received signal at the UAV from the BS can be expressed as

${y}=\sqrt{P_{b}G}{\bm{h}}^{\rm H}{\bm{\Theta}}{\bm{Hv}}u_{1}+{n}$ (2)

where $P_{b}$ is the source transmit power of BS,
${\rm\boldsymbol{v}}\in\mathbb{C}^{N\times 1}$ is the unit-magnitude precoding
vector, $u_{1}$ is the transmission signal with unit-power. ${n}$ is the
additive white Gaussian noise (AWGN) vector, following a distribution of
$\mathcal{CN}(0,\sigma_{0}^{2})$, and $\sigma_{0}^{2}$ is the noise power.
Then, the instantaneous signal-to-noise ratio (SNR) experienced by the UAV,
denoted by ${SNR}_{u}$, can be given by

${{SNR}_{u}=P_{b}G|\bm{h}^{\rm H}\bm{\Theta}\bm{Hv}|^{2}/\sigma_{0}^{2}}$ (3)

In the second transmission phase, the UAV decodes received URLLC packets from
the BS and re-encodes and forwards the packets to robots via the FDMA scheme.
Denote the UtG channel by $\left\\{g_{k}\right\\}$, $\forall
k\in{\mathcal{K}}$. Then the received signal of any robot in the considered
area can be expressed as

${y_{k}}=\sqrt{P_{u}}g_{k}{u_{2}}+{n_{k}}$ (4)

where $P_{u}$ is the transmit power of UAV, $u_{2}$ is the transmission signal
with unit-power, the noise ${n_{k}}$ is subject to a normal distribution with
zero mean and variance $\sigma_{k}^{2}$. Then, the instantaneous SNR
experienced by robot $k$, denoted by ${SNR}_{k}$, can be given by

${SNR}_{k}=P_{u}|g_{k}|^{2}/\sigma_{k}^{2}$ (5)

###  II-B BS-HAP Channel Model

In practice, HAP is usually deployed in an EM propagation environment without
rich scatters. It’s then reasonable to assume that the BS-HAP channel is
dominated by the LOS component.

The resulting BS-HAP channel matrix can be expressed as

${\bm{H}}=\alpha\bm{a}_{\rm RIS}(\upsilon_{{\rm AoD}}(\bm{p}_{H}))\bm{a}_{\rm
BS}^{\rm H}(\upsilon_{{\rm AoA}}(\bm{p}_{H}))$ (6)

where $\alpha$ represents the complex channel gain, capturing path loss and
additional attenuation caused by the EM propagation environment (e.g., rain
and water vapor). $\bm{a}_{\rm RIS}(.)\in\mathbb{C}^{M}$, $\bm{a}_{\rm
BS}(.)\in\mathbb{C}^{N}$ are the array steering vectors of BS and RIS-mounted
HAP, respectively, with

$\displaystyle\bm{a}_{\rm RIS}(.)=[1,e^{-j2\pi\bar{d}_{\rm RIS}{{\rm
sin}(.)}},...,e^{-j2\pi(M-1)\bar{d}_{\rm RIS}{{\rm sin}(.)}}]^{\rm T}$ (7)
$\displaystyle\bm{a}_{\rm BS}(.)=[1,e^{-j2\pi\bar{d}_{\rm BS}{{\rm
sin}(.)}},...,e^{-j2\pi(N-1)\bar{d}_{\rm BS}{{\rm sin}(.)}}]^{\rm T}$ (8)

where $\bar{d}_{\rm RIS}=d_{\rm RIS}/\lambda$, $\bar{d}_{\rm BS}=d_{\rm
BS}/\lambda$, $\lambda$ is wavelength. $\upsilon_{{\rm AoD}}(\bm{p}_{H})$ and
$\upsilon_{{\rm AoA}}(\bm{p}_{H})$ are the AoD and AoA of the BS-HAP link,
respectively. Aided by precise position information obtained by GPS or other
positioning systems, both AoA and AoD can be obtained.

Owing to the relative stability of the BS-HAP channel, we assume that its
channel gain can be maintained stable for a long time and continuous
estimation is not necessary. Therefore, the channel gain is considered to be
known in advance in this paper.

###  II-C HAP-UAV Channel Model

Unlike the BS-HAP channel, the EM propagation environment of the HAP-UAV
channel will be more complex. A UAV is usually deployed close to robots to
establish high-quality air-to-ground links; as a result, there is limited
physical scatters (e.g., trees, buildings, and mountains) around the UAV.
Accordingly, we assume that the HAP-UAV channel consists of limited
reflected/scattered paths. In addition, due to the high-frequency transmission
characteristic [24], the transmitted signal has a sound geometric
characteristic. To this end, we adopt a narrowband geometric block-fading
channel model in [25]

$\displaystyle\bm{h}=\sum\nolimits_{l=1}^{L}\beta_{l}\bm{a}_{\rm
RIS}(\omega_{l})$ (9)

where $\beta_{l}$ and $\omega_{l}$ represent the complex path gain coefficient
and AoD of the $l$-th path, respectively. We assume
$\beta_{l}\sim\mathcal{CN}(0,\sigma^{2})$, where $\sigma^{2}$ represents the
average power.

###  II-D UAV-to-Ground Channel Model

In the UtG channel model, we assume that the channel fading consists of
large–scale and small-scale channel fading. Specifically, the UtG channel can
be modeled as $\left\\{g_{k}=g_{k}^{L}g_{k}^{S}\right\\}$, where $g_{k}^{L}$
and $g_{k}^{S}$ represent large-scale and small-scale fading towards robot
$k$, respectively.

Large-scale channel statistics, including path loss and shadow fading, are
affected by the line-of-sight (LoS) probability of a UtG channel. We leverage
a widely adopted UtG 3-D channel model described in [26] to calculate large-
scale channel fading. In particular, a LoS component between the UAV and the
$k$-th robot exists with a specific probability, which is determined by the EM
propagation environment and locations of both the UAV and the GU.
Mathematically, the LoS probability, denoted by $P_{k}^{\rm
LoS}(\varphi_{k})$, $\forall k\in{\mathcal{K}}$, is given by [26]

$\displaystyle P_{k}^{\rm LoS}(\varphi_{k})=\frac{1}{1+\alpha_{1}{\rm
exp}(-\alpha_{2}(\varphi_{k}-\alpha_{1}))}$ (10)

where $\alpha_{1}$ and $\alpha_{2}$ are environment-dependent parameters,
$\varphi_{k}=\frac{180}{\pi}\times{\rm arctan}(\frac{h_{k}}{d_{k}})$ is the
elevation angle of the $k$-th robot towards the UAV, and
$d_{k}=\sqrt{(x_{U}-x_{G,k})^{2}+(y_{U}-y_{G,k})^{2}}$ is the horizontal
distance between the $k$-th robot and the UAV. Hence, the large-scale channel
gain can be expressed as

$\begin{array}[]{l}-20{\rm{lo}}{{\rm{g}}_{{\rm{10}}}}(g_{k}^{L})=20{\rm{lo}}{{\rm{g}}_{{\rm{10}}}}(4\pi{f_{\rm
UAV}}/c_{l})+20{\rm{lo}}{{\rm{g}}_{{\rm{10}}}}({d_{k}})+P_{k}^{{\rm{LoS}}}({\varphi_{k}}){\eta_{{\rm{LOS}}}}+P_{k}^{{\rm{NLoS}}}({\varphi_{k}}){\eta_{{\rm{NLOS}}}}\end{array}$
(11)

where $f_{\rm UAV}$ is the carrier frequency of UAV, ${c_{l}}$ is velocity of
light, $P_{k}^{\rm NLoS}(\varphi_{k})=1-P_{k}^{\rm LoS}(\varphi_{k})$ is the
probability of NLoS propagation. The environment-dependent parameters
$\eta_{\rm LOS}$ and $\eta_{\rm NLOS}$ are the excessive path-losses,
corresponding to the LoS and NLoS components, respectively.

The rayleigh fading model is utilized to model small-scale channel fading.
Further, the small-scale fading coefficient is assumed to be an independent
and identically distributed (i.i.d.) random variable with zero mean and unit
variance.

###  II-E Decoding Error Probability in a URLLC Regime

Different from conventional wireless communications, a URLLC regime has
stringent QoS requirements on transmission latency and reliability (e.g.,
DEP). To satisfy these requirements, the blocklength of URLLC packets should
be controlled within a finite length. However, the Shannon’s capacity formula
fails to build the relationship between these performance metrics and the
blocklength. Shannon’s capacity refers to the maximum achievable channel
coding rate (MAR) that can be obtained with arbitrary transmission latency,
which indicates that a signal blocklength should be sufficiently long. To
tackle this issue, we leverage tight MAR approximations [27, 28] of the BS-
HAP-UAV link (denoted by $R_{u}$) and the UAV-GU link (denoted by $R_{k}$),
which are extremely precise under the assumption that a blocklength is greater
than 50 symbols, i.e.,

$\displaystyle R_{*}={\rm
log_{2}}(1+{SNR}_{*})-\sqrt{\frac{1-(1+{SNR}_{*})^{-2}}{b_{*}}}\frac{Q^{-1}(\varepsilon_{*})}{\rm
ln2}$ (12)

where $*\in\left\\{u,k\right\\}$,
$Q(x)=\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}{\rm exp}(-t^{2}/2){\rm d}t$ is a
Gaussian Q-function, which is also a monotonically decreasing function of the
argument. $b_{u}$ and $b_{k}$ are the blocklengths allocated to the BS-HAP-UAV
link and the UAV-GU link, respectively. As the UAV adopts the DF strategy, the
BS-HAP-UAV link and the UAV-GU link can be configured with different
blocklengths. $\varepsilon_{u}$ and
$\left\\{\varepsilon_{k}\right\\}_{k=1}^{K}$ are the DEPs of the UAV and the
robots, respectively. Denote by $F_{k}$ the size of a URLLC packet transmitted
to the $k$-th robot from the UAV and $F_{B}$ the size of a URLLC packet sent
by the BS. According to the FDMA downlink transmission scheme, we have
$R_{u}=F_{B}/b_{u}$ and $R_{k}=F_{k}/b_{k}$ [29, 30, 31, 32]. The
corresponding DEP is then given by

$\begin{array}[]{l}{\varepsilon_{*}}=Q\left({\sqrt{\frac{{{b_{*}}}}{{1-{{(1+{{SN}}{{{R}}_{*}})}^{-2}}}}}\left({{\rm{lo}}{{\rm{g}}_{2}}\left({1+{{SN}}{{{R}}_{*}}}\right)-{R_{*}}}\right){\rm{ln}}2}\right)\buildrel\Delta\over{=}Q(f({b_{*}},{{SN}}{{{R}}_{*}}))\end{array}$
(13)

###  II-F Problem Formulation

Based on the above system model, this subsection aims to formulate a joint
transmission power, blocklength, and RIS phase shift optimization problem for
providing energy-efficient URLLC services for robots.

####  II-F1 URLLC and physical resource constraints

DEP and communication latency are crucial performance metrics of URLLC
services, which should be strictly constrained. Denote the DEPs of the UAV and
the $k$-th robot by $\varepsilon_{u}$ and $\varepsilon_{k}$, respectively. As
decoding errors, happening either at the UAV or at a GU, will result in
transmission failure, the overall DEP experienced by the $k$-th robot can be
written as

$\displaystyle\varepsilon=(1-\varepsilon_{u})\varepsilon_{k}+\varepsilon_{u}=\varepsilon_{k}+\varepsilon_{u}-\varepsilon_{u}\varepsilon_{k}<\varepsilon_{k}+\varepsilon_{u}$
(14)

Note that we omit the term $\varepsilon_{u}\varepsilon_{k}$ as it’s much
smaller than the other two terms. From (14), we can observe that reducing the
overall DEPs is equivalent to reducing the DEPs of UAV and robots,
respectively. What’s more, as a backbone link, it’s significant to provide an
extremely low and stable DEP over the BS-HAP-UAV link under any circumstances.
The statistical feature of the cascaded BS-HAP-UAV channel and the UtG channel
is also quite different. The cascaded BS-HAP-UAV channel is deterministic. The
UtG channel is random, and then an expected constraint on the DEP of a robot
should be enforced. Hence, we separately impose the DEP constraints of the UAV
and robots as follows,

$\varepsilon_{u}\leq\varepsilon^{\rm th}_{u}$ (15)
${\mathbb{E}}[\varepsilon_{k}]\leq\varepsilon^{\rm th}_{k},\forall
k\in{\mathcal{K}}$ (16)

where $\varepsilon^{\rm th}_{u}$ and $\varepsilon^{\rm th}_{k}$ are the
tolerable DEPs of the UAV and the $k$-th GU, respectively,
${\mathbb{E}}[\cdot]$ the expectation operator.

Next, we discuss the communication latency-related constraint. To guarantee
the ultra-low latency performance, the blocklength should be maintained at a
finite length [15, 30, 29]. Note that the UAV adopts the FDMA-DF mode, which
allows us to flexibly design different blocklengths over each communication
link to further boost the overall NSIN system performance. Define the minimum
and maximum blocklength sets allocated to the BS and the UAV as
$\left\\{B_{u}^{\min},B_{u}^{\max}\right\\}$ and
$\left\\{B^{\min},B^{\max}\right\\}$, respectively. Thus, the latency-related
constraints can be described as [30, 29]

$B_{u}^{\min}\leq b_{u}\leq B_{u}^{\max},\text{ }b_{u}\in\mathbb{Z^{+}}$ (17)
$B^{\min}\leq b_{k}\leq B^{\max},\text{ }b_{k}\in\mathbb{Z^{+}},\forall
k\in{\mathcal{K}}$ (18)

Finally, we discuss the constraints on BS and UAV transmit power. We
investigate the power control problem of BS and UAV to achieve energy-
efficient URLLC packet transmission. The configured transmit power is not
allowed to exceed a power budget. Denote the maximum transmit power of the BS
and the UAV by $P_{B}$ and $P_{U}$, respectively. Therefore, we have the
following transmit power constraints

$0<P_{b}\leq P_{B}$ (19) $0<P_{u}\leq P_{U}$ (20)

####  II-F2 Objective function design

Energy efficiency is a key performance evaluation indicator for UAV-assisted
cellular networks [33, 34, 35]. Based on the concepts of MAR, and BS and UAV
transmit power, we design an energy efficiency model to represent the
objective function. The energy efficiency model is defined as the ratio of the
effective throughput [29] to the system power consumption, which is given by

$\displaystyle\eta_{k}=\frac{R_{u}(1-\varepsilon_{u}/\varepsilon^{\rm
th}_{u})+R_{k}(1-\varepsilon_{k}/\varepsilon^{\rm
th}_{k})}{P_{b}/P_{B}+{P_{u}/P_{U}}},\forall k\in{\mathcal{K}}$ (21)

where the system power consumption denotes the sum of normalized BS and UAV
transmit power.

Remark 1: The value of DEP is usually extremely small, resulting in correct
decoding probabilities $1-\varepsilon_{u}$ and $1-\varepsilon_{k}$
approximately equals $1$. In this case, we must scale up the DEP value by
multiplying it by $\frac{1}{\varepsilon_{u}^{\rm th}}$ (or
$\frac{1}{\varepsilon_{k}^{\rm th}}$) to explicitly reflect the impact of
resource optimization on DEP. Additionally, as BS-HAP-UAV channel attenuation
is quite severe, the allocated BS transmit power will be very high to satisfy
stringent URLLC requirements. As a result, the BS transmit power will dominate
the total BS and UAV transmit power. Nevertheless, this case must be avoided
as the UAV is strictly energy-constrained and the UAV transmit power must be
carefully controlled. For this aim, we normalized the BS and UAV transmit
power by their maximum transmit power, respectively.

As different robots may experience diverse channel qualities, we aim at
maximizing the minimum system energy efficiency when guaranteeing the
stringent service requirements of all URLLC robots. We propose to jointly
optimize UAV and BS transmit power, blocklength, and RIS phase shift to
achieve the goal, under the constraints on DEPs, blocklengths, and UAV and BS
transmit power. Combined with the above analysis, the joint resource
optimization problem can be formulated as

$\displaystyle\max\limits_{P_{u},P_{b},b_{u},\\{b_{k}\\},\\{\phi_{n}\\}}\min\limits_{k\in\mathcal{K}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+\bar{p}_{u}}$
(22a) $\displaystyle{\rm s.t:}\text{
}Q\left({f\left({{b_{u}},{SNR}_{u}}\right)}\right)\leq\varepsilon_{u}^{{\rm{th}}}$
(22b)
$\displaystyle{{\mathbb{E}}[Q\left({f\left({{b_{k}},{SN}{{R}_{k}}}\right)}\right)]\leq\varepsilon_{k}^{{\rm{th}}},\forall
k\in{\mathcal{K}}}$ (22c) $\displaystyle\phi_{n}\in[0,2\pi),\forall
n\in\mathbb{N}$ (22d) $\displaystyle{\rm Constraints}\text{
}(\ref{eq_BS_length})-(\ref{eq_UAV_power_cons})$ (22e)

where ${\bar{p}}_{b}=P_{b}/P_{B}$, $\bar{p}_{u}=P_{u}/P_{U}$,
$\epsilon_{u}=\varepsilon_{u}/\varepsilon^{\rm th}_{u}$, and
$\epsilon_{k}=\varepsilon_{k}/\varepsilon^{\rm th}_{k}$. The solution of (22)
is highly challenging for the following reasons: 1) The HAP-UAV channel gain
is unknown and needs to be accurately estimated. However, there is a lack of
studies on sparse HAP-UAV channel modeling and corresponding accurate channel
estimation methods in the literature. 2) There are multiple EM propagation
paths in the HAP-UAV channel, and thus, it’s difficult to obtain the optimal
RIS phase shift such that the received signal strength of the UAV can be
maximized. 3) The expression of DEP is highly complex, which makes the problem
intractable and greatly hinders the theoretical analysis of the problem. 4)
Given the estimated HAP-UAV channel gain, (22) is still difficult to be
solved. This is because it contains integer and continuous variables and the
decision variables-related to HAP and UAV are intricately coupled, which
indicates that the problem is a mixed integer non-convex programming (MINCP)
problem.

In this paper, we develop a joint channel estimation, enhancement, and
resource optimization framework to solve the above challenging problem.In
particular, we first develop a probabilistic channel model to capture the
sparsity of the HAP-UAV channel and design a novel channel estimation approach
to estimate the HAP-UAV channel with high precision. Second, we design a RIS
phase shift optimization strategy that results in a suboptimal HAP-UAV channel
gain. Third, we derive the approximate expressions of DEPs to make the problem
analytically tractable. Fourth, we analyze the monotonicity of the objective
function of the problem. Based on the analysis results, we decompose the
problem into two-layered optimization problems, which reduces the difficulty
of the theoretical analysis of the problem.

##  III HAP-UAV CHANNEL ESTIMATION and ENHANCEMENT

In this section, we will introduce a message passing-based estimation approach
of the HAP-UAV channel. The RIS-user channel estimation is a hot and under-
studied topic. The recent works in [36, 37, 38] proposed to use compressed
sensing (CS) methods to estimate the RIS-user multi-path channel. However, the
authors in [36, 37, 38] did not discuss the crucial issue of angle estimation
or off-grid angle estimation. As a result, the accuracy of the channel
estimation needs to be improved. In this section, the HAP-UAV multi-path
channel with the off-grid angle consideration will be studied.

For multi-antenna communication systems, with an increasing number of
antennas, the spatial resolution in the angular domain will be extended
accordingly. As the HAP-UAV channel contains a finite number of propagation
paths, the channel will be very sparse in the angular domain. To effectively
utilize this sparsity property, we can transform the original channel model
into an angular domain representation. In this way, the channel estimation
problem is transformed into a sparse signal recovery problem. Generally, we
adopt a uniform sampling grid
${\bm{\widehat{\omega}}}=\left\\{\widehat{\omega}_{1},...,\widehat{\omega}_{N}\right\\}$,
which is obtained by discretizing the angular domain
$[-\frac{\pi}{2},\frac{\pi}{2}]$, to match the true AoDs denoted as
${\bm{\omega}}=\left\\{\omega_{1},...,\omega_{L}\right\\}$ ($L\ll N$). Under
the ideal conditions, all of the true AoDs can accurately match grid points.
However, the continuous distribution characteristic of the true angle
indicates that the grid points cannot effectively capture the true AoDs in
practice. Thus, the AoD estimation via an on-grid model usually leads to non-
negligible errors. To achieve high-precision estimation results, we introduce
off-grid offsets into the on-grid model, which is also called an off-grid
model. Specifically, given
$\omega_{l}\notin\left\\{\widehat{\omega}_{1},...,\widehat{\omega}_{N}\right\\}$
and its nearest sampling grid point
$\widehat{\omega}_{nl}=\mathop{\rm\min}\limits_{n_{l}\in\left\\{1,...,N\right\\}}\left\\{|\omega_{l}-\widehat{\omega}_{n_{l}}|\right\\}$,
we can define the off-grid offset $\Delta\phi_{nl}$ as

$\Delta\phi_{nl}=\begin{cases}{\omega}_{l}-\widehat{\omega}_{nl},&l\in{\mathcal{L}}=\\{1,2,\ldots,L\\}\\\
0,&{\rm otherwise}\end{cases}$ (23)

The case of $\Delta\phi_{nl}=0$ indicates that only active paths will be
assigned with a non-zero offset (If we don’t know the number of true paths in
advance, we can set a threshold of path strength to filter non-active paths).
Then the corresponding off-grid steering matrix is
${{\bm{A}(\widehat{\bm{\omega}}+\Delta\bm{\phi})}}=[\bm{a}(\widehat{\omega}_{1}+\Delta\phi_{1}),...,\bm{a}(\widehat{\omega}_{N}+\Delta\phi_{N})]$,
and the sparse channel vector in the complex angular domain is
$\bm{x}\in\mathbb{C}^{N}$ with $L$ non-zero elements. Finally, the off-grid
sparse angular domain representation of the downlink channel
$\displaystyle\bm{h}\in\mathbb{C}^{N}$ can be expressed as

$\bm{h}={\bm{A}(\widehat{\bm{\omega}}+\Delta\bm{\phi})\bm{x}}\triangleq{\bm{A}(\Delta\bm{\phi})\bm{x}}$
(24)

It’s noteworthy that there must be some off-grid offsets making (24) hold.
Though it’s quite challenging to find the optimal off-grid offsets, the gap
between the sampling grid points and the true AoDs will be significantly
reduced by iteratively adjusting the off-grid offsets. Therefore, the off-grid
model always outweighs the on-grid model.

###  III-A Probability Model of Channel Vector

Recall that the HAP-UAV channel has limited paths in the angular domain. To
fully exploit this feature, the probabilistic model of
${\bm{x}}=[x_{1},...,x_{N}]^{\rm T}$ can be modeled by two independent random
hidden vectors, i.e., _hidden channel support vector_ and _hidden channel
value vector_. The _hidden channel support vector_
${\bm{d}}=[d_{1},...,d_{M}]^{\rm T}\in\left\\{0,1\right\\}^{N}$ represents
whether the paths are active or not. Specifically, if $d_{n}=1$, then there is
an active path around the $n$-th AoD direction
$\widehat{\omega}_{N}+\Delta\delta_{N}$. We define $p(d_{n}=1)=\lambda$,
indicating the channel sparsity. The complex-valued _hidden channel value
vector_ ${\bm{q}}=[q_{1},...,q_{N}]^{\rm T}$ represents the complex gains of
paths. $q_{n}\sim\mathcal{CN}(\zeta,\rho)$ follows an i.i.d. complex white
Gaussian distribution with a mean $\zeta$ and variance $\rho$. Note that
$q_{n}$ does not determine whether the $n$-th path from the corresponding AoD
direction is active or not. The channel vector can then be expressed as

${\bm{x}}=[d_{1}q_{1},...,d_{N}q_{N}]^{\rm T}$ (25)

Accordingly, the joint probabilistic HAP-UAV channel model with a prior
distribution can be formulated as

$p({\bm{x}})=p({\bm{x}},{\bm{d}},{\bm{q}})=p({\bm{x}}|{\bm{d}},{\bm{q}})p({\bm{d}})p({\bm{q}})$
(26)

where $p({\bm{x}}|{\bm{d}},{\bm{q}})$ is the joint conditional prior. From
(26), we know that $p(x_{n}|d_{n},q_{n})=1$ conditioned on $x_{n}=d_{n}q_{n}$;
otherwise, $p(x_{n}|d_{n},q_{n})=0$. Thus, $p({\bm{x}}|{\bm{d}},{\bm{q}})$ can
be expressed as

$\displaystyle
p({\bm{x}}|{\bm{d}},{\bm{q}})=\prod_{n=1}^{N}p(x_{n}|d_{n},q_{n})=\prod_{n=1}^{N}\delta(x_{n}-d_{n}q_{n})$
(27)

where $\delta(.)$ is the Dirac delta function.

###  III-B Problem Formulation of Off-Grid Sparse Channel Estimation

Using the angular domain channel representation, we can rewrite (2) as a
standard CS model

${{\bm{y}}={\bm{F}}(\Delta\bm{\phi}){\bm{x}}+\bm{n}_{e}}$ (28)

where ${\bm{F}}(\Delta\bm{\phi})={\bm{U_{0}}^{\rm H}}{\bm{\bar{H}}}^{\rm
H}{\bm{\Theta}}^{\rm H}{\bm{A}(\Delta\bm{\phi})}\in\mathbb{C}^{P\times N}$ is
the measurement matrix, ${\bm{\bar{H}}}$ is the normalized BS-HAP channel
matrix, ${\bm{U_{0}}}\in\mathbb{C}^{M\times P}$ is a training pilot matrix
transmitted from the BS, $\bm{n}_{e}\sim\mathcal{CN}(0,\sigma_{e}^{2}\bm{I})$
is the complex AWGN. For a CS problem, the choice of measurement matrix
directly affects the performance of the recovery approach. In the proposed
model, we design pilot $\bm{U}_{0}$ based on the partial discrete Fourier
transform (DFT) random permutation (pDFT-RP) measurement matrix, whose
effectiveness is confirmed in [40]. The following Lemma presents how to
generate the pilot matrix ${\bm{U_{0}}}$.

######  Lemma 1.

For any off-grid offset $\Delta\phi$, the pilot matrix to be transmitted by
the BS can be expressed as

${{\bm{U}_{0}^{\rm H}}={\bm{SDR}}({{{\bm{\bar{H}}}^{\rm H}}{{\bm{\Theta}}^{\rm
H}}\bm{A}(\bm{0})})^{-1}}$ (29)

###### Proof.

Please refer to Appendix A in the technical report [39]. ∎

We can then estimate the off-grid sparse channel with (28). Besides, we
leverage the minimum mean square error (MMSE) criterion to evaluate the
estimation performance. The estimated complex gain is
$\hat{x}_{n}=E(x_{n}|{\bm{y}},\Delta\bm{\phi})$, which will result in the
MMSE. The expectation is taken over the marginal posterior
$p(x_{n}|{\bm{y}},\Delta\bm{\phi})$ given measurements ${\bm{y}}$ and off-grid
offsets $\Delta\bm{\phi}$. To exploit the inherent statistical structure of
the HAP-UAV channel in terms of channel support and value vectors, the
marginal posterior distribution can be expressed by Bayes’ rule as

$\begin{array}[]{l}\displaystyle
p(x_{n}|{\bm{y}},\Delta\bm{\phi})\propto\sum_{\bm{d}}\int_{{\bm{x}}_{-n},{\bm{q}}}p({\bm{x,d,q}|{{\bm{y}},}\Delta\bm{\phi}})=\sum\limits_{\bm{d}}\int_{{\bm{x}}_{-n},{\bm{q}}}p({\bm{d}})p({\bm{q}})\prod_{n=1}^{N}p(x_{n}|d_{n},q_{n})\prod_{p=1}^{P}p(y_{p}|{\bm{x}},\Delta\bm{\phi})\end{array}$
(30)

where $\propto$ indicates proportionality up to a constant scale factor,
${\bm{x}}_{-m}$ means a vector $\bm{x}$ excluding the element $x_{m}$.
According to (28), the likelihood function
$p(y_{p}|{\bm{x}},\Delta\bm{\phi})=\mathcal{CN}(y_{p};{\bm{f}}_{p}{\bm{x}},\sigma_{e}^{2})$,
where ${\bm{f}}_{p}$ is the $p$-th row of measurement matrix $\bm{F}$. The
optimal off-grid offsets $\Delta\bm{\phi}$ can be obtained by a maximum
likelihood (ML) approach

$\Delta\bm{\phi}^{*}={\rm arg}\max\limits_{\Delta\bm{\phi}}{\rm
ln}p({\bm{y}},\Delta\bm{\phi})={\rm arg}\max\limits_{\Delta\bm{\phi}}{\rm
ln}\int_{\bm{x}}p({\bm{x}},{\bm{y}},\Delta\bm{\phi}){\rm d}{\bm{x}}$ (31)

After obtaining the estimated offsets $\Delta\bm{\phi}^{*}$, we can then
calculate the corresponding marginal posterior distribution (30). It’s the
production of many distribution functions, each of which is determined by only
a subset of variables. Such decomposition structures can be calculated by a
graphical model, i.e., a factor graph, which can be an undirected bipartite
graph used to represent the connection between random variables (called
variable nodes) and related probability density functions (pdfs) (called
factor nodes). Unfortunately, in (30), the generated factor graph contains
densely connected structure due to the measurement matrix, which will greatly
hinder the exact calculation of the posterior distribution using the factor
graph. To handle this challenge, we propose a Recurrent-OAMP approach that can
decouple the linear and nonlinear observations in the following subsection.

###  III-C Recurrent-OAMP Approach

The Recurrent-OAMP approach contains two iteratively working modules, as shown
in Fig. 3: Module A is a linear MMSE (LMMSE) estimator, inputting observations
and messages from Module B; Module B is an MMSE estimator, inputting the
inherent sparse prior and messages from Module A. Both modules will work
iteratively until convergence.

Figure 2: Two recurrent modules of the R-OAMP approach.

Figure 3: A factor graph of hidden channel support and value vectors.

####  III-C1  Module A: LMMSE Estimator

In Module A, since only linear observations are considered, we can consider
that $\bm{x}$ has a Gaussian distribution with mean ${\bm{x}}^{pri}_{A}$ and
variance $v_{A}^{pri}$ passed from module B without a sparse prior. Given
received signals $\bm{y}$ and the prior distribution
$\mathcal{CN}(\bm{x};{\bm{x}}^{pri}_{A},v_{A}^{pri}\bm{I})$, the LMMSE
estimation ${\bm{x}}^{post}_{A}$ and its mean squared error (MSE)
${v_{A}^{post}}$ are respectively given by [40]

$\displaystyle{\bm{x}}^{post}_{A}={\bm{x}}^{pri}_{A}+\frac{v_{A}^{pri}}{v_{A}^{pri}+\sigma_{e}^{2}}{\bm{F}(\Delta\bm{\phi})}^{\rm
H}({\bm{y}}-{\bm{F}(\Delta\bm{\phi})}{\bm{x}}^{pri}_{A})$ (32)
$\displaystyle{v_{A}^{post}}={v_{A}^{pri}}-\frac{P}{N}.\frac{(v_{A}^{pri})^{2}}{v_{A}^{pri}+\sigma_{e}^{2}}$
(33)

As $\mathcal{CN}({\bm{x}};{\bm{x}}^{post}_{A},v_{A}^{post}\bm{I})$ and
$\mathcal{CN}({\bm{x}};{\bm{x}}^{pri}_{A},v_{A}^{pri}\bm{I})$ are independent
of each other, the extrinsic message is also subject ot a Gaussian
distribution and given by

$\mathcal{CN}({\bm{x}};{\bm{x}}^{post}_{A},v_{A}^{post}\bm{I})\propto\mathcal{CN}({\bm{x}};{\bm{x}}^{pri}_{A},v_{A}^{pri}\bm{I})\mathcal{CN}({\bm{x}};{\bm{x}}^{ext}_{A},v_{A}^{ext}\bm{I})$
(34)

Hence, the extrinsic mean and variance can be computed by [40]

$\displaystyle{\bm{x}}^{pri}_{B}={\bm{x}}^{ext}_{A}={v}^{ext}_{A}(\frac{{\bm{x}}^{post}_{A}}{{v}^{post}_{A}}-\frac{{\bm{x}}^{pri}_{A}}{{v}^{pri}_{A}})$
(35)
$\displaystyle{v}^{pri}_{B}={v}^{ext}_{A}=(\frac{1}{{v}^{post}_{A}}-\frac{1}{{v}^{pri}_{A}})^{-1}$
(36)

where ${\bm{x}}^{pri}_{B}$ and ${v}^{pri}_{B}$ are the input mean and variance
of Module B, which are equal to ${\bm{x}}^{ext}_{A}$ and ${v}^{ext}_{A}$,
respectively.

####  III-C2 Module B: MMSE Estimator

In Module B, the MMSE estimator is designed with proposed sparsity channel
priors. However, an MMSE estimator in the OAMP approach executes estimation
based on i.i.d. priors, which thus cannot be directly applied. Therefore, we
need to redesign the MMSE estimator by fully exploiting the sparsity structure
in the HAP-UAV channel. For the message passing-based sparse signal recovery
approach, a basic assumption is that the extrinsic message
${\bm{x}}^{pri}_{B}$ is achieved by an approximated additive white Gaussian
noise (AWGN) observation of ${\bm{x}}$, mathematically,

$\displaystyle{\bm{x}}^{pri}_{B}={\bm{x}}+{\bm{z}}$ (37)

where ${\bm{z}}\sim\mathcal{CN}(0,{v}^{pri}_{B}\bm{I})$ is independent of
${\bm{x}}$. Extensive experiments have proved that the AWGN observation
assumption is valid for OAMP [41]. Additionally, if constructing a factor
graph using (37) rather than the joint distribution
$p({\bm{x,d,q}{,{\bm{y}}}})$, the per iteration computational complexity of a
message passing approach running on the factor graph can be significantly
reduced to $O(N)$. Therefore, based on (37), the factor graph of the joint
distribution $p({\bm{x,d,q}{,{\boldsymbol{x}}^{pri}_{B}}})$ denoted by
$\mathcal{J}$ can be shown in Fig. 3, where the functional form of each factor
node is listed as below: $g_{n}(x_{n},x^{pri}_{B,n})\sim
p(x^{pri}_{B,n}|x_{n})=\mathcal{CN}({x_{n}};x^{pri}_{B,n},v_{B}^{pri})$,
$f_{n}(x_{n},d_{n},q_{n})\sim p(x_{n}|d_{n},q_{n})=\delta(x_{n}-d_{n}q_{n})$,
$u_{n}(d_{n})\sim p(d_{n})={(1-\lambda)^{1-d_{n}}(\lambda)^{d_{n}}}$,
$e_{n}(q_{n})\sim p(q_{n})=\mathcal{CN}({q_{n}};\zeta,\rho)$. Our goal is to
calculate the marginal posterior distribution
$\left\\{p({x_{n}}|{\bm{x}}^{pri}_{B})\right\\}_{n=1}^{N}$ by message passing
on the graph $\mathcal{J}$. As the graph $\mathcal{J}$ does not contain loops,
we can directly explore sum-product message passing rules without additional
approximation. The details of the message scheduling are as below.

Step 1: According to the sum-product algorithm, the message from factor node
$g_{n}$ to variable node $x_{n}$ is

$\eta_{g_{n}\rightarrow
x_{n}}(x_{n})=p(x^{pri}_{B,n}|x_{n})=\mathcal{CN}({x_{n}};x^{pri}_{B,n},v_{B}^{pri})$
(38)

Step 2: Given prior information, compute the message from variable node
$q_{n}$ to factor node $f_{n}$, i.e.,

$\eta_{q_{n}\rightarrow f_{n}}(q_{n})=\mathcal{CN}({q_{n}};\zeta,\rho)$ (39)

Step 3: According to sum-product rules,the message from factor node $f_{n}$ to
variable node $x_{n}$ is

$\begin{array}[]{l}{\eta_{{f_{n}}\to{x_{n}}}}({x_{n}})=\frac{{\sum\limits_{{d_{n}}}{\int\limits_{q_{n}}{{\eta_{{q_{n}}\to{f_{n}}}}}}(q_{n}){\eta_{{d_{n}}\to{f_{n}}}}(d_{n})\delta(x_{n}-{d_{n}}q_{n}){\rm
d}q_{n}}}{{\sum\limits_{{d_{n}}}{\int\limits_{x_{n}}{\int\limits_{q_{n}}{{\eta_{{q_{n}}\to{f_{n}}}}}}}(q_{n}){\eta_{{d_{n}}\to{f_{n}}}}(d_{n})\delta(x_{n}-{d_{n}}q_{n}){\rm
d}q_{n}{\rm
d}x_{n}}}=(1-\lambda)\delta(x_{n})+\lambda\mathcal{CN}(\delta(x_{n});\zeta,\rho)\end{array}$
(40)

It’s noteworthy that (40) is a Bernoulli-Gaussian distribution. With just
three parameters $\left\\{\lambda,\zeta,\rho\right\\}$, we can control the
distribution of $x_{n}$. In practice, even though we can’t predict these
parameters, we can estimate their values by some methods (e.g., expectation-
maximization).

Step 4: Then, the posterior distribution can be given by

$p({x_{n}}|{\bm{x}}^{pri}_{B})\propto\eta_{f_{n}\rightarrow
x_{n}}(x_{n})\eta_{g_{n}\rightarrow x_{n}}(x_{n})$ (41)

According to the principle of MMSE, the posterior mean and variance can be
respectively computed by

$\displaystyle x^{post}_{B,n}={\rm
E}({x_{n}}|{\bm{x}}^{pri}_{B})=\int_{x_{n}}{x_{n}}p({x_{n}}|{\bm{x}}^{pri}_{B})$
(42) $\begin{array}[]{l}\displaystyle
v^{post}_{B}=\frac{1}{N}\sum_{n=1}^{N}{\rm
Var}({x_{n}}|{\bm{x}}^{pri}_{B})=\frac{1}{N}\sum_{n=1}^{N}\int_{x_{n}}|{x_{n}}-{\rm
E}({x_{n}}|{\bm{x}}^{pri}_{B})|^{2}p({x_{n}}|{\bm{x}}^{pri}_{B})\end{array}$
(43)

After updating the posterior distribution of ${\bm{x}}$, the extrinsic update
mean and covariance passed to Module A can be expressed as

$\displaystyle{\bm{x}}^{pri}_{A}={\bm{x}}^{ext}_{B}={v}^{ext}_{B}(\frac{{\bm{x}}^{post}_{B}}{{v}^{post}_{B}}-\frac{{\bm{x}}^{pri}_{B}}{{v}^{pri}_{B}})$
(44)
$\displaystyle{v}^{pri}_{A}={v}^{ext}_{B}=(\frac{1}{{v}^{post}_{B}}-\frac{1}{{v}^{pri}_{B}})^{-1}$
(45)

When the posterior distribution is obtained by message passing in Module B,
the extrinsic message will be sent back to module A. This message-passing
pattern is similar to how a recurrent network works, so the recurrent-OAMP
(R-OAMP) approach is named. R-OAMP executes message passing between Module A
and Module B until convergence.

Once the estimated posterior distribution
$\hat{p}({\bm{x}|{\bm{y},\Delta\bm{\phi}}})=\prod_{n=1}^{N}p(x_{n}|{\bm{y}},\Delta\bm{\phi})$
is obtained, we can update the off-grid offsets $\Delta\bm{\phi}$ using ML.
Unfortunately, we can’t obtain the closed-form expression of the log-
likelihood function ${\rm ln}p({\bm{y}},\Delta\bm{\phi})$. As a result, the
gradient descent method can’t be directly utilized to gain the optimal
solution $\Delta\bm{\phi}^{*}$. To address this challenging issue, we utilize
the in-exact majorization-minimization (MM) algorithm in [42] to approximate
the optimal solution. As a general form of the expectation-maximization
algorithm, the key idea of the in-exact MM algorithm is to iteratively
generate a progressively refined (log-likelihood) function (a lower bound of
${\rm ln}p({\bm{y}},\Delta\bm{\phi})$) and approximate the optimal solution by
continuously improving the lower bound. Although the expectation-maximization
algorithm can’t guarantee the theoretical optimality, it can guarantee to
converge to a stationary point [43]. Specifically, let the surrogate function
of ${\rm ln}p({\bm{y}},\Delta\bm{\phi})$ be
$r(\Delta\bm{\phi};\Delta\bm{\dot{\phi}})$, which indicates that this
progressively refined function is constructed at a given point
$\Delta\bm{\dot{\phi}}$. As the lower bound,
$r(\Delta\bm{\phi};\Delta\bm{\dot{\phi}})$ should satisfy the following three
conditions:

$r(\Delta\bm{\phi};\Delta\bm{\dot{\phi}})\leq{\rm
ln}p({\bm{y}},\Delta\bm{\phi}),\forall\Delta\bm{\phi}$ (46)
$r(\Delta\bm{\dot{\phi}};\Delta\bm{\dot{\phi}})={\rm
ln}p({\bm{y}},\Delta\bm{\dot{\phi}})$ (47) $\displaystyle\frac{\partial
r(\Delta\bm{\phi};\Delta\bm{\dot{\phi}})}{\partial\Delta\bm{\phi}}|_{\Delta\bm{\phi}=\Delta\bm{\dot{\phi}}}=\frac{\partial{\rm
ln}p(\bm{y},\Delta\bm{\phi})}{\partial\Delta\bm{\phi}}|_{\Delta\bm{\phi}=\Delta\bm{\dot{\phi}}}$
(48)

It’s not hard to verify that the following constructed progressively refined
function satisfies all the above conditions

$\displaystyle r(\Delta\bm{\phi};\Delta\bm{\dot{\phi}})=\int
p({\bm{x}}|{\bm{y}},\Delta\bm{\dot{\phi}}){\rm
ln}\frac{p({\bm{x}},{\bm{y}},\Delta\bm{\phi})}{p({\bm{x}}|{\bm{y}},\Delta\bm{\dot{\phi}})}d{{\bm{x}}}$
(49)

With the constructed function, we can iteratively update $\Delta\bm{\phi}$ by

$\Delta\bm{\phi}^{i+1}={\rm
arg}\max\limits_{\Delta\bm{\phi}}r(\Delta\bm{\phi};\Delta\bm{\phi}^{i})$ (50)

where $i$ represents the $i$-th iteration and $\Delta\bm{\phi}^{i}$ stands for
the corresponding value of $\Delta\bm{\phi}$.

Nevertheless, it’s quite difficult to obtain the optimal
$\Delta\bm{\phi}^{i+1}$ as (50) is non-convex. To solve this problem
effectively, we explore a gradient ascend method, and $\Delta\bm{\phi}^{i+1}$
can be updated by

$\displaystyle\Delta\bm{\phi}^{i+1}=\Delta\bm{\phi}^{i}+\alpha_{i}\frac{\partial
r(\Delta\bm{\phi};\Delta\bm{\phi}^{i})}{\partial\Delta\bm{\phi}}|_{\Delta\bm{\phi}=\Delta\bm{\phi}^{i}}$
(51)

where $\alpha_{i}$ is the update stepsize. Besides, to obtain high-precision
estimation results, we utilize the backtracking linear search method, which
can avoid oscillating results or slow convergence by adaptively adjusting the
iteration stepsize.

Please refer to Appendix B in the technical report [39]. for the detailed
update expression of off-grid offsets $\Delta\bm{\phi}^{i+1}$ .

Algorithm 1 Recurrent-OAMP

1: Input: Received signal (or observations) $\bm{y}$, measurement matrix
${\bm{F}}(0)$, noise variance $\sigma^{2}_{e}$, $I_{c}^{\max}=50$.

2: Initialize: ${\bm{x}}^{pri}_{A}=\bm{0}$, ${v}^{pri}_{A}=\lambda\rho$,
$\Delta\bm{\phi}=\bm{0}$

3: while not converge do

4: Module A (given $\Delta\bm{\phi}=\Delta\bm{\phi}^{i}$):

5: Update ${\rm\boldsymbol{x}}^{post}_{A}$ and ${v}^{post}_{A}$ using (32) and
(33).

6: Update ${\bm{x}}^{pri}_{B}={\bm{x}}^{ext}_{A}$ and
${v}^{pri}_{B}={v}^{ext}_{A}$ using (35) and (36).

7: Module B:

8: Message passing over the path $g_{n}\rightarrow x_{n}$ using (38) with
outputs of Module A.

9: Message passing over the path $q_{n}\rightarrow f_{n}$ using (39).

10: Message passing over the path $f_{n}\rightarrow x_{n}$ using (40).

11: Calculate the posterior distribution $f_{n}\rightarrow x_{n}$ using (41),
and update ${\bm{x}}^{post}_{B}$ and ${v}^{post}_{B}$ using (42) and (43).

12: Update ${\bm{x}}^{pri}_{A}={\bm{x}}^{ext}_{B}$ and
${v}^{pri}_{A}={v}^{ext}_{B}$ using (44) and (45).

13: Repeat Module A and Module B until convergence.

14: Output
$\hat{p}({\bm{x}|{\bm{y},{\Delta\bm{\phi}}^{i}}})=\mathcal{CN}({x_{n}};{\bm{x}}^{post}_{B},v_{B}^{post}\bm{I})$

15: Off-grid Estimation:

16: Update off-grid offsets $\Delta\bm{\phi}^{i+1}$ using (84).

17: Update $i=i+1$.

18: end while

19: Output: $\hat{\bm{x}}={\bm{x}}_{B}^{post}$ and
$\hat{\bm{h}}={\bm{A}(\Delta\bm{\phi}^{i+1})\hat{\bm{x}}}$.

###  III-D Channel Enhancement via RIS Phase Alignment

In this subsection, we discuss how to enhance the HAP-UAV channel by aligning
RIS phases. The optimization of RIS phase shift is crucial for the performance
improvement of RIS-integrated communication systems. By aligning the phase
shift ${\bm{\Theta}}$, signals sent by the BS and then reflected by the RIS
can be coherently added to the UAV, which thus significantly increase the HAP-
UAV channel gain and maximize the received SNR by the UAV ${SNR}_{u}$.

Given the estimated HAP-UAV channel vector $\bm{h}$, we can then formulate the
RIS phase shift optimization problem as
${\bm{\Theta}}^{\star}=\mathop{\arg\max}\limits_{\bm{\Theta}}\frac{P_{b}G}{\sigma_{0}^{2}}|\bm{h}^{\rm
H}{\bm{\Theta}}{\bm{H}}\bm{{v}}|^{2}$, where $\sigma_{0}^{2}$ is the noise
power. To achieve the optimal ${\bm{\Theta}}$, we first need to determine the
precoding vector $\bm{v}$. The following Lemma presents a strategy for
designing the precoding vector.

######  Lemma 2.

For any given BS-HAP-UAV propagation channel and RIS phase shift, the BS needs
to adopt a maximum ratio transmission (MRT) strategy to maximize ${SNR}_{u}$,
that is, we have

$\displaystyle\bm{v}=\frac{\bm{a}_{BS}(\upsilon_{{\rm
AoD}})}{||\bm{a}_{BS}(\upsilon_{{\rm AoD}})||_{2}}$ (52)

###### Proof.

Please refer to Appendix C in the technical report [39]. ∎

Using (52), ${SNR}_{u}$ can be rewritten as

$\begin{array}[]{l}SN{R_{u}}=\frac{{{P_{b}}G}}{{{\sigma_{0}^{2}}}}\left|{{{\bm{h}}^{\rm{H}}}{\bm{\Theta}}{\bm{H}}\frac{{{{\bm{a}}_{{\rm{BS}}}}({\upsilon_{{\rm{AoD}}}})}}{{||{{\bm{a}}_{{\rm{BS}}}}({\upsilon_{{\rm{AoD}}}})|{|_{2}}}}}\right|=\frac{{\alpha{P_{b}}G}}{{{\sigma_{0}^{2}}}}{\left|{\sum\limits_{n=1}^{N}{\sum\limits_{l=1}^{L}{{\beta_{l}}}}{e^{j({\theta_{n}}+2\pi(n-1){{\bar{d}}_{{\rm{RIS}}}}({\rm{sin}}({\omega_{l}})-{\rm{sin}}({\upsilon_{{\rm{AoA}}}}))}}}\right|^{2}}\end{array}$
(53)

where $\beta_{l}$ is the $l$-th non-zero element of the estimated gain
$\hat{\bm{x}}$.

From (53), we can find that ${SNR}_{u}$ is independent of AoD $\upsilon_{{\rm
AoD}}$ of BS-HAP link. However, as there are multiple reflected/scattered
paths in the HAP-UAV EM propagation environment, it’s difficult the optimize
$\theta_{n}$ ($\forall n$) such that the signal received power of the UAV can
be maximized. Hence, we attempt to obtain the suboptimal $\theta_{n}$
($\forall n$). To this end, we rearrange (53) as

$\displaystyle{SNR}_{u}=\frac{\alpha
P_{b}G}{\sigma_{0}^{2}}\left|\sum_{l=1}^{L}\beta_{l}r_{l}\right|^{2}$ (54)

where $r_{l}\triangleq|\sum_{n=1}^{N}e^{j(\theta_{n}+2\pi(n-1)d_{\rm RIS}({\rm
sin}(\omega_{l})-{\rm sin}(\upsilon_{{\rm AoA}}))}|^{2}$.

Therefore, maximizing (54) is equivalent to maximizing
$\\{r_{l}\\}_{l\in{\mathcal{L}}}$. Denote $l_{0}$ as the phase alignment path
of RIS and $\psi_{n}^{*}(l_{0})$ the corresponding phase shift. Then, we can
coherently reflect signals via path $l_{0}$ by setting

$\psi_{n}^{*}(l_{0})=\theta-2\pi(n-1){\bar{d}}_{\rm RIS}({\rm
sin}(\omega_{l_{0}})-{\rm sin}(\upsilon_{{\rm AoA}}))$ (55)

where $\theta\in[0,2\pi)$ is an arbitrary phase shift. By substituting
$\theta_{n}$ for $\psi_{n}^{*}(l_{0})$ and performing specific mathematical
transformations, we can then rewrite $r_{l}$ as

$\displaystyle r_{l}=g(\Delta\psi_{l}(l_{0}))\triangleq\left|\frac{{\rm
sin}(\pi Nd_{\rm RIS}\Delta\psi_{l}(l_{0}))}{{\rm sin}(\pi d_{\rm
RIS}\Delta\psi_{l}(l_{0}))}\right|^{2}$ (56)

where $\Delta\psi_{l}(l_{0})={\rm sin}(\omega_{l})-{\rm sin}(\omega_{l_{0}})$.
To maximize (56), we need to minimize $|\Delta\psi_{l}(l_{0})|$. Considering
the low-latency requirement of URLLC services, we perform the first-order
Taylor expansion on $|\Delta\psi_{l}(l_{0})|$ to realize the rapid
configuration of RIS phase shift, i.e.,

$|\Delta\psi_{l}(l_{0})|\approx|{\rm
cos}(\omega_{l})||\omega_{l}-\omega_{l_{0}}|$ (57)

As a result, we can reformulate the RIS Phase shift optimization problem as

$\displaystyle\min\limits_{\omega_{l_{0}},\zeta_{l}}\sum_{l=1}^{L}\left|{\beta_{l}}{{\rm
cos}(\omega_{l})}\right|\zeta_{l}$ (58a) $\displaystyle{\rm s.t:}\text{
}\omega_{l_{0}}-\zeta_{l}\leq\omega_{l},\forall l\in{\mathcal{L}}$ (58b)
$\displaystyle\omega_{l}\leq\omega_{l_{0}}+\zeta_{l},\forall
l\in{\mathcal{L}}$ (58c)

(58) is a linear programming problem and can then be efficiently solved by
optimization tools such as MOSEK.

##  IV RESOURCE OPTIMIZATION and ALGORITHM DESIGN

We next discuss how to solve the original optimization problem with the
cascaded BS-HAP-UAV channel gain. As mentioned above, the complicated DEP
expression greatly hinders the theoretical analysis of the optimization
problem. Therefore, we first attempt to approximate the DEP.

Figure 4: The function ${\varepsilon_{u}}$ in comparison with its two types of
approximations.

###  IV-A DEP Approximation

To obtain an analysis-friend expression of DEP ${\varepsilon_{u}}$, we improve
the linear approximation method in [44]. The key idea of this method is to
conduct the first-order Taylor series expansion at the point
${SNR}_{u}=\gamma_{u}$ to achieve an asymptotical expression. Different from
the approximation in [44], we extend the linear approximation region and
propose to allocate redundant resources to achieve a lower DEP, as shown in
Fig. 4. Specifically, we approximate $\varepsilon_{u}$ as

$\begin{array}[]{l}{\varepsilon_{u}}\approx\Omega({SNR}_{u})=\left\\{{\begin{array}[]{*{20}{l}}{1,}&{SNR_{u}\leq{{SNR}}_{\rm
low}}\\\ {0,}&{SNR_{u}\geq{{SNR}}_{\rm up}}\\\
{\frac{1}{2}-\frac{1}{{2\chi_{u}}}({{SN}}{{{R}}_{u}}-{\gamma_{u}}).}&{\rm
otherwise}\end{array}}\right.\end{array}$ (59)

where $\displaystyle\gamma_{u}=2^{R_{u}}-1$, ${SNR}_{\rm
low}={{\gamma_{u}}-\frac{1}{{{\chi_{u}}}}}$, ${SNR}_{\rm
up}={{\gamma_{u}}+\frac{1}{{{\chi_{u}}}}}$, and
$\displaystyle\chi_{u}=\sqrt{\frac{b_{u}}{2\pi}}(2^{2R_{u}}-1)^{-1/2}$.
Similarly, $\Omega({SNR}_{k})$ represents the Taylor series expansion based
approximation of $Q(f(b_{k},{SNR}_{k}))$.

Recall that the UtG channel is a random channel. However, the UtG channel
fading should be regarded as a constant during a transmission interval. Thus,
we need to analyze the expected DEP, denoted by
${\overline{\varepsilon}_{k}}$, experienced by any robot $k$, which can be
calculated by

$\begin{array}[]{l}{\overline{\varepsilon}_{k}}={\mathbb{E}}[Q(f({b_{k}},SN{R_{k}}))]\approx\int_{0}^{\infty}\Omega({SNR}_{k})f_{{SNR}_{k}}(x){\rm
d}x=\chi_{k}\int_{\gamma_{k}-\frac{1}{{\chi_{k}}}}^{\gamma_{k}+\frac{1}{{\chi_{k}}}}F_{{SNR}_{k}}(x){\rm
d}x\end{array}$ (60)

where $f_{{SNR}_{k}}(x)$ and $F_{{SNR}_{k}}(x)$ are the pdf and the cdf of
${SNR}_{k}$, respectively.

As the Rayleigh small-scale fading is considered, ${SNR}_{k}$ follows an
exponential distribution. Thus, $F_{{SNR}_{k}}(x)$ can be computed by

$F_{{SNR}_{k}}(x)=\begin{cases}1-e^{-\frac{x}{\overline{{SNR}}_{k}}},&x\geq
0\\\ 0,&x<0\end{cases}$ (61)

where $\overline{{SNR}}_{k}=P_{u}|g_{k}^{L}|^{2}/\sigma_{k}^{2}$ is the
average value of ${SNR}_{k}$. By substituting (61) into (60), for any robot
$k\in{\mathcal{K}}$, we can calculate its expected DEP by

$\begin{array}[]{l}{\overline{\varepsilon}_{k}}\approx{\chi_{k}}\int_{{\gamma_{k}}-\frac{1}{{{\chi_{k}}}}}^{{\gamma_{k}}+\frac{1}{{{\chi_{k}}}}}{(1-{e^{-x/{{\overline{SNR}}_{k}}}})}{\rm
d}x{\mathop{\approx}\limits^{(a)}\frac{2\gamma_{k}}{{\overline{SNR}}_{k}}}\end{array}$
(62)

where (a) follows from the approximation
$1-e^{-x/\overline{{SNR}}_{k}}\approx\frac{x}{\overline{{SNR}}_{k}}$ under an
assumption of high SNR [45]. This assumption is reasonable for URLLC
applications because robots need to achieve high SNR such that the stringent
QoS requirements can be satisfied.

###  IV-B Original Problem Decomposition

After obtaining the estimated HAP-UAV channel gain, the aligned RIS phase, and
the approximated DEPs, the original problem (22) becomes analytically
tractable. Nevertheless, we observe that there are integer and continuous
variables in (22) and the decision variables are complexly coupled. (22) can
be confirmed to be a MINCP problem, which is challenging. Generally, resorting
to an iterative optimization scheme is a feasible way of solving this
challenging problem. Yet, the developed algorithm would contain multiple
nested loops if we directly adopt the iterative optimization scheme, which
increases the computational complexity of the algorithm.

Fortunately, given a BS blocklength $b$, we observe that the optimal $P_{b}$
does not vary with the UAV transmit power $P_{u}$. This is because there is a
narrow feasible region
$\hat{\mathcal{D}}_{p}\subseteq{\mathcal{D}_{p}}=\left\\{{{P_{b}}\left|{Q\left({f\left({{b_{u}},SN{R_{u}}}\right)}\right)\leq\varepsilon_{u}^{{\rm{th}}},0<{P_{b}}\leq{P_{B}}}\right.}\right\\}$,
in which $\Omega(SNR_{u})$ strictly monotonically decreases with $P_{b}$.
Given the blocklength $b$, the objective function (22a) then reduces to a
monotonically increasing function over $P_{b}$ in the specific feasible region
$\hat{\mathcal{D}}_{p}$.

Therefore, we can apply horizontal decomposition to (22) that creates a
problem of two-layered structure. The $1^{\rm st}$ layer is _BS-Layer
Optimization_ , and the $2^{\rm nd}$ layer is _UAV-Layer Optimization_. More
importantly, we need not to conduct iterative optimization between these two
layers.

####  IV-B1 BS-Layer Optimization

Given any decision variables-related to the UAV-Layer Optimization, we can
formulate the BS-Layer Optimization problem as

$\displaystyle\max\limits_{P_{b}}\text{ }\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+{\bar{p}}_{u}}$
(63a) $\displaystyle{\rm s.t:}\text{ }\Omega(SNR_{u})\leq\varepsilon_{u}^{\rm
th}$ (63b) $\displaystyle{\rm Constraints\text{
}(\ref{eq_BS_length}),(\ref{eq_BS_power})}$ (63c)

where $\min\left\\{x_{k}\right\\}$ denotes the minimum $x_{k}$, $\forall
k\in{{\mathcal{K}}}$.

We observe from the expressions of $R_{u}$ and $\Omega(SNR_{u})$ that the
decision variables $P_{b}$ and $b_{u}$ in (63) are intricately coupled, which
challenges the solution to (63). To tackle this issue, we propose to decompose
(63) into two subproblems, named _BS transmit power control_ subproblem and
_BS blocklength optimization_ subproblem.

1-a) BS transmit power control: Given $b_{u}$, ${P}_{u}$ and $\\{b_{k}\\}$,
the BS power control subproblem can be reduced to

$\displaystyle\max\limits_{P_{b}}\text{ }\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+{\bar{p}}_{u}}\text{
}{\rm s.t:}\text{ }{\rm Constraints\text{
}(\ref{eq_BS_power}),(\ref{eq:BS_layer_opt}b)}$ (64a)

Recall that $\Omega({SNR}_{u})$ is a piecewise function. Thus, we need to
discuss the power control subproblem in three different cases closely related
to the value of ${SNR}_{u}$. We can directly exclude the case of
${SNR}_{u}\leq{\gamma_{u}}-\frac{1}{{\chi_{u}}}$ as it will result in an
unacceptable DEP. Then, we discuss the remaining two cases, respectively.

Case I: ($SNR_{u}\geq SNR_{\rm up}$) In this case, we have
$\varepsilon_{u}=0$, and the BS power control subproblem can be reformulated
as

$\displaystyle\max\limits_{P_{b}}\text{ }\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+{\bar{p}}_{u}}$
(65a) $\displaystyle{\rm s.t:}\text{ }\displaystyle{\bar{p}}_{b}\geq{SNR_{\rm
up}}/({P_{B}\Delta B})$ (65b) $\displaystyle{\rm Constraint\text{
}(\ref{eq_BS_power}).}$ (65c)

where $\Delta B=G|\bm{h}^{\rm H}\bm{\Theta}\bm{Hv}|^{2}/\sigma_{0}^{2}$. As
the numerator of (65a) is not related to ${\bar{p}}_{u}$, the objective
function of (65) monotonically decreases with ${P}_{b}$. Thus, given $b_{k}$,
the optimized BS transmit power (denoted by $P_{b}(b_{k})$) can be computed by

$\displaystyle P_{b}(b_{u})={\rm\min}\left\\{P_{B},{SNR_{\rm up}}/{\Delta
B}\right\\}$ (66)

Case II: ($SNR_{\rm low}<SNR_{u}<SNR_{\rm up}$) In this case, We can
reformulate (64) as

$\displaystyle\max\limits_{P_{b}}\text{ }\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+{\bar{p}}_{u}}$
(67a) $\displaystyle{\rm s.t:}\text{ }\displaystyle{SNR_{\rm
low}}/({P_{B}\Delta B})\leq{\bar{p}}_{b}\leq{SNR_{\rm up}}/({P_{B}\Delta B})$
(67b) $\displaystyle{\rm Constraint\text{ }(\ref{eq_BS_power}).}$ (67c)

Generally, the monotonicity of (67a) is hard to be determined. Fortunately, we
observe that the feasible region ${\mathcal{P}}_{b}$ of (67) is rather narrow,
in which the numerator of (67a) varies greatly, while the denominator remains
almost unchanged. Therefore, the monotonicity of (67a) is dominated by its
numerator. The numerator is then a monotonically increasing function in
${\mathcal{P}}_{b}$, and so is (67a). It further indicates that the two
discussed cases have the same solution that is given in (66).

1-b) BS blocklength optimization problem: Given $P_{b}(b_{u})$, $P_{u}$, and
$\\{b_{k}\\}$, we can formulate the BS blocklength optimization problem as

$\displaystyle\max\limits_{b_{u}}\text{
}\min\limits_{k\in{\mathcal{K}}}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{P_{b}(b_{B})/P_{B}+\overline{p}_{u}}\text{
}{\rm s.t:}\text{ }{\rm Constraints\text{
}(\ref{eq_BS_length}),(\ref{eq:BS_layer_opt}b).}$ (68a)

Note that the obtained $P_{b}(b_{u})$ is a function of $b_{u}$. After
substituting $P_{b}(b_{u})$ into (68), we can obtain an objective function
regarding variable $b_{u}$. However, the obtained objective function is non-
concave even given a relaxed variable $b_{u}$. Fortunately, considering that
the feasible region of (68) is small, it will be a good choice to achieve its
optimal solution by exhaustively searching.

Remark 2: Observe that the optimized $P_{b}(b_{u})$ is a function $b_{u}$.
After obtaining $b_{u}$ using an exhaustive search method, we can directly
obtain the optimized $P_{b}(b_{u})$ by (66). It indicates that we can solve
the _BS-Layer Optimization_ problem in just two iterations.

####  IV-B2 UAV-Layer Optimization

Given any decision variables related to the BS-Layer Optimization, we can
formulate the UAV-Layer Optimization problem as

$\displaystyle\max\limits_{P_{u},\\{b_{k}\\}}\text{
}\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+{\bar{p}}_{u}}$
(69a) $\displaystyle{\rm s.t:}\text{
}{{2\gamma_{k}}}/{{\overline{SNR}}_{k}}\leq\varepsilon_{k}^{\rm th},\forall
k\in{\mathcal{K}}$ (69b) $\displaystyle{\rm Constraints\text{
}(\ref{eq_UE_length}),(\ref{eq_UAV_power_cons})}$ (69c)

Similarly, as (69) is a mixed-integer optimization problem, we decompose it
into two iteratively optimized subproblems, i.e., _UAV transmit power control_
subproblem and _UAV blocklength optimization_ subproblem, such that (69) can
be effectively solved.

2-a) UAV transmit power control: Assuming that the values of the blocklength
of $b_{k}$, $\forall k$, ${P}_{b}$, and $b_{u}$ are given, we can then
formulate the UAV transmit power control subproblem as

$\displaystyle\max\limits_{P_{u}}\text{ }\eta_{\rm
EE}=\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+{\bar{p}}_{u}}\text{
}{\rm s.t:}\text{ }{\rm Constraints\text{
}(\ref{eq_UAV_power_cons}),(\ref{eq:UAV_layer_opt}b).}$ (70a)

The objective function of (70) is a complex fraction, which makes it hard to
be directly optimized. To address this issue, we define the optimal $\eta_{\rm
EE}^{\star}$ as

$\displaystyle\eta_{\rm EE}^{\star}=\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}({\bar{p}}_{u}^{\star})(1-{\epsilon}_{k}({\bar{p}}_{u}^{\star}))}{{\bar{p}}_{b}+{\bar{p}}_{u}^{\star}}$
(71)

where ${\bar{p}}_{u}^{\star}=P_{u}^{\star}/{P_{U}}$ is the optimal solution of
(70). The following Lemma presents how to calculate $\eta_{\rm EE}^{\star}$.

######  Lemma 3.

$\eta_{\rm EE}^{\star}$ can be achieved if and only if

$\begin{array}[]{l}\mathop{\max}\limits_{{P_{u}}}\mathop{\min}\limits_{k\in{\mathcal{K}}}\text{
}R_{u}(1-\epsilon_{u})+{R_{k}}(1-\epsilon{{}_{k}})-{\eta_{\rm
EE}^{\star}}\left({{{\bar{p}}_{b}}+{{\bar{p}}_{u}}}\right)\\\
=\mathop{\min}\limits_{k\in{\mathcal{K}}}\text{
}R_{u}(1-\epsilon_{u})+{R_{k}}({{\bar{p}}_{u}^{\star}})(1-\epsilon{{}_{k}}({{\bar{p}}_{u}^{\star}}))-{\eta_{\rm
EE}^{\star}}\left({{{\bar{p}}_{b}}+{{\bar{p}}_{u}^{\star}}}\right)\\\
=0\end{array}$ (72)

###### Proof.

Please refer to Appendix D in the technical report [39]. ∎

According to Lemma 3, we can transform the expression of the objective
function of (70). As a result, (70) can be equivalently transformed into

$\displaystyle\max\limits_{P_{u}}\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle R_{u}(1-\epsilon_{u})+R_{k}(1-{\epsilon}_{k})-\eta_{\rm
EE}({{\bar{p}}_{b}+{\bar{p}}_{u}})$ (73a) $\displaystyle{\rm s.t:}\text{
}\displaystyle{\bar{p}}_{u}\geq{{2\gamma_{k}}\sigma_{k}^{2}}/({|g_{k}^{L}|^{2}\varepsilon^{\rm
th}_{k}P_{U}}),\forall k\in{\mathcal{K}}$ (73b) $\displaystyle{\rm
Constraint\text{ }(\ref{eq_UAV_power_cons}).}$ (73c)

Note that (73) is a typical max-min optimization problem. Then, we introduce
an auxiliary variable $y$ and equivalently transform (73) into

$\displaystyle\max\limits_{P_{u},y}\text{ }y$ (74a) $\displaystyle{\rm
s.t:}\text{
}\displaystyle{{2\gamma_{k}}\sigma_{k}^{2}}/({|g_{k}^{L}|^{2}{\bar{p}}_{u}P_{U}})\leq
1+({R_{u}(1-\epsilon_{u})-y-\eta_{\rm
EE}({{\bar{p}}_{b}+{\bar{p}}_{u}})})/{R_{k}},\forall k\in{\mathcal{K}}$ (74b)
$\displaystyle{\rm Constraints\text{
}(\ref{eq_UAV_power_cons}),(\ref{eq:reformu_UAV_power_control}b).}$ (74c)

It can be confirmed that (74) is a rotated quadratic cone programming problem
that can be efficiently optimized by some commercial optimization tools such
as MOSEK. Therefore, we can conclude the steps of optimizing (70) in the
Algorithm 2.

Algorithm 2 UAV Transmit Power Control

1: Initialization: Initialize $\eta_{\rm EE}^{{(0)}}$ and $y^{(0)}$. Let error
tolerance parameter $\Delta_{r}=1{\rm e}-3$, maximum iterations $r_{\rm
max}=50$, and $r=0$. Input $\\{b_{k}\\}$ and the optimized $b_{u}$ and
$P_{b}$.

2: repeat

3: Given $\eta_{\rm EE}^{{(r)}}$, solve (74) to obtain the optimized UAV
transmit power ${P}_{u}^{(r+1)}$ using MOSEK.

4: Given ${P}_{u}^{(r+1)}$ and $\eta_{\rm EE}^{{(r)}}$, calculate $y^{(r+1)}$
by  ${y^{(r+1)}}=\mathop{\min}\limits_{k\in{\cal K}}\text{
}{R_{u}}(1-\epsilon{{}_{u}})+{R_{k}}(1-{{2{\gamma_{k}}\sigma_{k}^{2}}}/({{|{g_{k}}{|^{2}}\bar{p}_{u}^{(r+1)}{P_{U}}}}))-\eta_{{\rm{EE}}}^{(r)}\left({{{\bar{p}}_{b}}+\bar{p}_{u}^{(r+1)}}\right)$
(75)

5: Given $y^{(r+1)}$ and ${P}_{u}^{(r+1)}$, update $\eta_{\rm EE}^{{(r+1)}}$
by  $\displaystyle\eta_{\rm
EE}^{{(r+1)}}={\min\limits_{k\in{\mathcal{K}}}\left\\{R_{u}(1-\epsilon_{u})+R_{k}(1-{{2\gamma_{k}}\sigma_{k}^{2}}/({|g_{k}|^{2}{\bar{p}}_{u}^{(r+1)}P_{U}})))\right\\}}/({{\bar{p}}_{b}+{\bar{p}}_{u}^{(r+1)}})$
(76)

6: until The convergence criterion $|y^{(r+1)}-y^{(r)}|\leq\Delta_{r}$ is met
or reach the maximum number of iterations $r_{\rm max}$.

2-b) UAV blocklength optimization Given the UAV transmit power $P_{u}$, we can
formulate the UAV blocklength optimization subproblem as

$\displaystyle\max_{\\{b_{k}\\}}\min\limits_{k\in{\mathcal{K}}}\text{
}R_{k}(1-\epsilon_{k})\text{ }{\rm s.t:}\text{ }\rm{Constraints\text{
}(\ref{eq_UE_length}),(\ref{eq:UAV_layer_opt}b)}$ (77a)

Note that (77) is a non-linear integer programming problem. The computational
complexity of solving it using some traditional optimization methods is
extremely high. Nevertheless, the feasible region of (77) is small. Then, we
can resort to the exhaustive search to achieve its optimal solution quickly.
Moreover, we can search $\\{b_{k}\\}_{k=1}^{K}$ in parallel.

###  IV-C Algorithm Design

Finally, according to the above analysis and derivation, we can summarize the
main steps of the joint RIS Phase shift, Transmit Power, and Blocklength
optimization (PTPB) algorithm for mitigating the formulated energy efficiency
optimization problem in Algorithm 3.

Algorithm 3 Joint RIS Phase shift, Transmit Power, and Blocklength
optimization, PTPB

1: Initialization: Run the initialization steps in R-OAMP and Algorithm 2. Let
$c_{\rm max}=100$, $r=0$. $I_{\rm max}=10$.

2: Channel estimation and RIS phase shift optimization:

3: Call the channel estimation approach (R-OAMP) to estimate the HAP-UAV
channel gain $\bm{h}$.

4: Compute $\bm{v}$ by (52) and optimize (58) using MOSEK to obtain
$\omega_{l_{0}}$. Let $\theta_{n}=\psi_{n}^{*}(l_{0})$.

5: Given the obtained $\bm{h}$, $\bm{v}$, and $\bm{\Theta}$, compute $SNR_{u}$
by (3).

6: BS-Layer optimization:

7: repeat

8: Given $b_{u}^{(i)}$ compute $P_{b}(b_{u}^{(i)})^{(i+1)}$ by (66).

9: Given $P_{b}(b_{u}^{(i)})^{(i+1)}$, and any non-zero $P_{U}$ and
$\\{b_{k}\\}$, exhaustively search $b_{u}^{(i+1)}$ such that (68a) is
maximized, subject to (68b). Update $i=i+1$.

10: until Converge or reach the maximum number of iterations $I_{\rm max}$.

11: UAV-Layer optimization:

12: repeat

13: Call Algorithm 2 to obtain $P_{u}^{(c+1)}$.

14: for each $k\in\\{1,2,\ldots,K\\}$ in parallel do

15: Given $P_{u}^{(c+1)}$, solve the following problem using an exhaustive
search method  $\displaystyle\max_{{b_{k}}}\text{ }R_{k}(1-\epsilon_{k})$
(78a) $\displaystyle{\rm s.t:}\text{
}{{2\gamma_{k}}}/{{\overline{SNR}}_{k}}\leq\varepsilon_{k}^{\rm th}$ (78b)
$\displaystyle B^{\min}\leq b_{k}\leq B^{\max},\text{ }b_{k}\in\mathbb{Z^{+}}$
(78c)

16: end for

17: Update $c=c+1$.

18: until Converge or $c=c_{\rm max}$.

The computational complexity of the proposed algorithm consists of three
contributors: 1) _Channel estimation_ : R-OAMP includes a module A, a module
B, and a parameter updating module. Module A is an LMMSE estimator with a
complexity of $O(PN)$. Module B is an MMSE estimator where a sum-product
approach is leveraged, which thus has a low computational complexity of
$O(N)$. The complexity of updating offsets is $O(PN^{2})$. Considering that
the above modules should be iteratively executed for no more than
$I_{c}^{\max}$ times. The computational complexity of R-OAMP is
$O(f_{1})=O(I_{c}^{\max}(PN+N+PN^{2}))$ in the worst-case. However, empirical
evidence shows that R-OAMP usually converges in 30 iterations. 2) _RIS phase
shift optimization_ : The complexity of solving (58) using an interior-point
method dominates the complexity of this contributor and is $O((1+L)^{3.5})$
[46]. 3) _Joint power and blocklength optimization_ : This joint optimization
model includes the BS-Layer optimization and the UAV-Layer optimization. As
the iterative method in BS-Layer optimization will converge in two iterations,
the complexity of BS-Layer optimization is dominated by the exhaustive search
approach and is
$O(f_{2})=O((B_{u}^{\max}-B_{u}^{\min})S_{1}+\log_{2}(B_{u}^{\max}-B_{u}^{\min}))$,
where $S_{1}$ is the complexity of computing (68a). For UAV-Layer
optimization, it needs to iteratively optimize (70) and (77) for no more than
$c_{\max}$ iterations. The complexities of optimizing UAV transmit power and
blocklength are $O(r_{\max}2^{3.5})$ and
$O((B^{\max}-B^{\min})S_{2}+\log_{2}(B^{\max}-B^{\min}))$ per iteration, where
$S_{2}$ is the complexity of computing (77a). Thus, the complexity of UAV-
Layer optimization is
$O(f_{3})=O(c_{\max}(r_{\max}2^{3.5}+(B^{\max}-B^{\min})S_{2}+\log_{2}(B^{\max}-B^{\min})))$
in the worst-case. Then, the total computational complexity of is
$O(f_{1}+f_{2}+(1+L)^{3.5}+f_{3})$ in the worst-case.

##  V Simulation results

###  V-A Comparison Algorithms and Simulation Parameters

In this section, we evaluate the performance of the proposed channel
estimation approach, RIS phase shift optimization strategy, and resource
optimization algorithm, respectively.

For this aim, the following benchmark algorithms are respectively compared:
OMP [47]: The OMP approach is an improved version of the well-known matching
pursuit approach. SBL [42]: The sparse Bayesian learning (SBL) approach
recovers sparse channel vectors based on a sparse Bayesian learning theory. SP
[48]: The subspace pursuit (SP) approach identifies the signal support based
on the maximum correlation criterion and is robust to measurement noises.
Random-Phase: In this algorithm, each element of the HAP-RIS is assigned with
a random phase shift. Zero-Phase: In this algorithm, all elements of the HAP-
RIS are assigned with a zero phase shift. Exhaustive-Phase: It exhaustively
searches for the AoD of the optimal phase alignment path ${\omega_{l_{0}}}$ in
the interval $[0,\pi/2]$. MTP: The difference between the maximum transmit
power-based (MTP) algorithm and the proposed algorithm lies in that MTP adopts
the maximum BS transmit power and the maximum UAV transmit power. MBL: The
difference between the maximum blocklength-based (MBL) algorithm and the
proposed algorithm lies in that MBL is configured with the maximum BS
blocklength and the maximum UAV blocklength.

We consider a square area $\mathcal{G}$ of size $500\times 500$ m2, the
location of the center point of which is $[80\ {\rm km},0\ {\rm km}]$. The
robots are randomly distributed in $\mathcal{G}$. The carrier frequency of BS
$f_{\rm BS}=6$ GHz, and the carrier frequency of UAV $f_{\rm UAV}=2$ GHz. The
BS-HAP channel attenuation and the variance $\sigma^{2}$ are generated by the
propagation model in[49]. The channel realization consists of $L=8$
propagation paths to ensure the sparsity characteristic, and the AoDs
$\\{\omega_{l}\\}$ are concentrated in an angular interval of ${\pi/12}$
centered on the LoS direction between the HAP and the UAV. The inter-element
spacing of the RIS is assumed to be half-wavelength. More system parameters
are listed as below: $M=32$, $N=128$, ${\bm{p}}_{H}=[18,1,0]^{\rm T}$ km,
$K=10$, ${\bm{p}}_{U}=[0.05,80,0]^{\rm T}$ km, $F_{B}=80$ bits,
$B^{\min}_{u}=100$ bits, $B^{\min}=100$ bits, $B^{\max}_{u}=1000$ bits,
$B^{\max}=1000$ bits, $P_{B}=120$ W, $P_{U}=0.5$ W, $\varepsilon^{\rm
th}_{u}=0.00005$, $\varepsilon^{\rm th}_{k}=0.00005$, $\alpha_{1}=12.08$,
$\alpha_{2}=0.11$, $\eta_{\rm LOS}=1.6$ dB, $\eta_{\rm NLOS}=23$ dB, $G=4$ dB,
$\sigma_{0}^{2}=-134$ dBm, $\sigma^{2}_{k}=-143$ dBm, $F_{k}=80$ bits [26,
30].

###  V-B Performance Evaluation

To comprehensively evaluate the performance of the proposed algorithm, we
design the simulation of validating the accuracy of the proposed channel
estimation approach, the simulation of verifying the effectiveness of the
designed RIS phase shift optimization strategy, and the simulation of
testifying the superiority of the proposed resource allocation algorithm,
respectively. All comparison approaches and algorithms in the simulations are
executed 200 times, and the final simulation result is the averaged one.

####  V-B1 Results of Channel Estimation

In this simulation, we compare the proposed R-OAMP approach with other channel
estimation approaches, including OMP, SBL, and SP. The parameters
$\left\\{{\lambda,\varsigma,\rho}\right\\}$ are automatically updated by
adopting an expectation-maximization algorithm during each iteration [50]. The
metric, normalized mean square error (NMSE), is selected to evaluate the
performance of channel estimation approaches, which is defined as

$NMSE(\hat{\bm{h}})\buildrel\Delta\over{=}{{{{\|{\hat{\bm{h}}-{\bm{h}}}\|}^{2}}}}/{{{{\left\|{\bm{h}}\right\|}^{2}}}}$
(79)

where ${\hat{\bm{h}}}$ is the estimation of ${\bm{h}}$.

We plot the tendency of estimation accuracy of all compared channel estimation
approaches under different $SNR_{u}\in\\{0,4,8,12,16,20\\}$, as shown in Fig.
5. From this figure, we can observe that: 1) For all channel estimation
approaches, their estimation performance will be improved with an increasing
$SNR_{u}$. 2) Compared with other estimation approaches, the estimation
accuracy is improved by at least $1.46$ dB using R-OAMP. 3) When there are
more RIS reflecting elements, which indicates that the HAP-UAV channel becomes
sparser, the estimation accuracy of the proposed approach is further improved.

(a) $N=128$

(b) $N=256$

Figure 5: NMSE versus received SNR of UAV ${\rm SNR}_{u}$, $P=48$.

The estimation performance of all compared channel estimation approaches
versus the number of pilot sequences $P$ is illustrated in Fig. 6. We can
obtain the following observations from this figure: 1) Except for OMP, the
obtained NMSE of the rest estimation approaches decreases sharply with an
increasing number of pilot sequences. For example, the obtained NMSE by R-OAMP
decreases by $15.24$ dB when $P$ increases from $30$ to $70$. 2) R-OAMP can
always achieve the smallest NMSE under different $P$ and $N$. It demonstrates
that R-OAMP can achieve more accurate channel estimation results with lower
pilot overhead by exploiting the sparsity prior of the HAP-UAV channel in the
angular domain.

(a) $N=128$

(b) $N=256$

Figure 6: NMSE versus the number of pilot sequences $P$, ${\rm SNR}_{u}=15\rm
dB$.

####  V-B2 Results of RIS Phase Shift Optimization

In this simulation, we compare the proposed algorithm with two benchmarks,
i.e., Random-Phase and Zero-Phase, to verify its effectiveness regarding RIS
phase shift optimization.

The tendency of the cascaded BS-HAP-UAV channel gain versus the number of RIS
passive reflecting elements $N\in\\{96,128,160,192,224,256\\}$ is plotted in
Fig. 8. We can observe from this figure that: 1) PTPB outweighs the benchmarks
and improves the BS-HAP-UAV channel gain by at least $20.36$ dB. 2) PTPB
achieves a cascaded channel gain that is close to the channel gain obtained by
the Exhaustive-Phase algorithm, and their difference is less than $3$ dB. 3)
The obtained channel gain by PTPB monotonically increases with the number of
RIS elements, while the benchmarks fail to unlock the significant advantage of
RIS in terms of addressing the channel fading issue.

Figure 7: BS-HAP-UAV channel gain versus the number of RIS elements $N$.

Figure 8: Achieved energy efficiency versus the number of RIS elements $N$.

Besides, we adopt the evaluation metric, energy efficiency computed by (22a),
to further verify the effectiveness of the designed RIS phase shift
optimization strategy. We plot the obtained energy efficiency of the
comparison algorithms under different numbers of RIS elements and reliability
constraints in Fig. 8. By referring to Fig. 8, the achieved BS-HAP-UAV channel
gains by both Random-Phase and Zero-Phase are extremely small. The results in
Fig. 8 illustrate that the proposed algorithm can achieve higher energy
efficiency and satisfy the reliability requirement with a lower transmit power
budget. From the results in Figs. 8 and 8, we can draw the conclusion that
PTPB can effectively deal with the phase alignment issue in a RIS-integrated
multipath EM propagation environment.

####  V-B3 Results of Transmit Power and Blocklength Optimization

To verify the effectiveness of the proposed algorithm regarding joint transmit
power and blocklength optimization, we compare it with two benchmark
algorithms, i.e., MTP and MBL. To this end, we evaluate the impact of the
number of the RIS elements $N$, the side length of $\mathcal{G}$, and the UAV
transmit power budget $P_{U}$, respectively.

(a) $\varepsilon^{\rm th}_{u}=\varepsilon^{\rm th}_{k}=0.00005$

(b) $\varepsilon^{\rm th}_{u}=\varepsilon^{\rm th}_{k}=0.000005$

Figure 9: Achieved energy efficiency versus the number of RIS elements $N$.

Fig. 9 shows the achieved energy efficiency of all comparison algorithms over
the number of RIS elements. From this figure, we can obtain the following
observations: 1) PTPB outperforms the other two benchmarks under diverse $N$,
$\varepsilon^{\rm th}_{u}$, and $\varepsilon^{\rm th}_{k}$. The achieved
energy efficiency by all comparison algorithms increases with an increasing
number of RIS elements. This is due to the improvement in the BS-HAP-UAV
channel quality, and then fewer resources (including transmit power and
blocklength) are required to satisfy URLLC requirements. 2) When a more
stringent reliability requirement, i.e., $\varepsilon^{\rm
th}_{u}=\varepsilon^{\rm th}_{k}=0.000005$, is enforced, the obtained energy
efficiency of all comparison algorithms will decrease. Obviously, more
resources will be consumed to satisfy the more stringent reliability
requirement. 3) Under a more stringent reliability constraint, the energy
efficiency obtained by MBL does not increase effectively with improved channel
conditions. In this case, the transmit power cannot be significantly reduced
owing to the strictly limited blocklength.

(a) $\varepsilon^{\rm th}_{u}=\varepsilon^{\rm th}_{k}=0.00005$

(b) $\varepsilon^{\rm th}_{u}=\varepsilon^{\rm th}_{k}=0.000005$

Figure 10: Achieved energy efficiency versus the side length of $\mathcal{G}$.

Fig. 10 plots the tendency of the obtained energy efficiency of all comparison
algorithms versus the side length of the coverage area $\mathcal{G}$. From
this figure, we can observe that: 1) PTPB obtains the maximum energy
efficiency under different side lengths of $\mathcal{G}$. 2) The achieved
energy efficiency by PTPB decreases as the area of the coverage area
increases. Enlarging the area of the coverage area indicates that more robots
will be distributed at the edge of the UAV cell, which thus consumes more
resources and decreases the achieved energy efficiency.

(a) $\varepsilon^{\rm th}_{u}=\varepsilon^{\rm th}_{k}=0.00005$

(b) $\varepsilon^{\rm th}_{u}=\varepsilon^{\rm th}_{k}=0.000005$

Figure 11: Achieved energy efficiency versus the transmit power budget of UAV.

Finally, we plot the tendency of the obtained energy efficiency of all
comparison algorithms over the transmit power budget of UAV
$P_{U}\in\\{400,500,600,700,800,900,1000\\}$ mW in Fig. 11. In this
simulation, we set the maximum blocklength $B_{u}^{max}=B^{max}=1000$ bits.
The following observations can be obtained from this figure: 1) PTPB achieves
the highest energy efficiency. For instance, under a tighter reliability
constraint, the system energy efficiency can be improved by at least $34.14\%$
by executing the proposed PTPB algorithm. 2) The performance of MBL degrades
more compared to MTP when a more stringent reliability constraint is imposed.
The reasons are as follows: compared to the case of $\varepsilon^{\rm
th}_{u}=\varepsilon^{\rm th}_{k}=0.00005$, MBL needs to allocate higher
transmit power to satisfy the tighter reliability requirement in the case of
$\varepsilon^{\rm th}_{u}=\varepsilon^{\rm th}_{k}=0.000005$. Owing to the
limited feasible region of blocklength and the adoption of the maximum
transmit power, the obtained energy efficiency by MTP will not be
significantly reduced when a more stringent reliability requirement is
imposed.

Summarily, from the above simulation results, we can conclude that joint RIS
phase shift, transmit power, and blocklength optimization can significantly
improve the system energy efficiency.

##  VI Conclusion

This paper proposed to deploy a RIS-integrated NSIN to provide energy-
efficient URLLC services for remote robots, which leveraged the advantages of
RIS and power control. For this aim, we formulated the energy-efficient URLLC
service provision problem as a resource optimization problem aiming at
maximizing the effective throughput and minimizing the system power
consumption, subject to URLLC and physical resource constraints. This problem
was difficult to be mitigated due to the unknown channel model and intractable
analysis. To handle these challenges, we developed a novel channel estimation
approach, derived an analysis-friend expression of DEP, and decomposed the
problem by analyzing its monotonicity. Based on the above results, we proposed
a joint RIS phase shift, transmit power, and blocklength optimization
algorithm, the computational complexity of which was also analyzed. Simulation
results verified that the developed channel estimation approach could
accurately estimate HAP-UAV channels and the proposed resource optimization
algorithm was more energy-efficient than benchmarks. This paper investigated a
stationary scenario, and extending it to a mobility scenario (e.g., the UAV is
moving) will be an interesting topic worthy of study in the near future. In
this case, channel tracking and RIS phase shift optimization approaches
adapted to moving UAVs need to be designed.

###  -A Proof of Lemma 1

In practice, we don’t know the off-grid offsets in advance. Nevertheless, if
we deploy a great number of antenna elements, the measurement matrix
${\bm{F}}(\Delta\bm{\phi})\approx{\bm{F}}(\boldsymbol{0})$. Therefore, without
loss of generality, we set the corresponding off-grid offsets to zero when
generating the pilot matrix. Particularly, we then have

${\bm{F}}(\Delta\bm{\phi})\approx{\bm{F}}({\bm{0}})=\bm{SDR}$ (80)

where ${\bm{S}}=\left\\{0,1\right\\}^{P\times N}$ is a selection matrix
generated by randomly selecting and reordering $P$ rows of an ${N\times N}$
identity matrix, ${\bm{D}}\in\mathbb{C}^{N\times N}$ is the DFT matrix,
${\bm{R}}\in\mathbb{C}^{N\times N}$ is a random permutation matrix generated
by randomly reordering an ${N\times N}$ identity matrix. Note that the
measurement matrix comprises the pilot matrix, BS-HAP channel matrix, and RIS
phase shift matrix. We can compute the channel state information (CSI) of the
BS-HAP channel in advance and can impose $\theta_{n}=-2\pi(n-1){\bar{d}}_{\rm
RIS}({\sin}({\arctan(\frac{||[x_{H},y_{H}]^{\rm T}-[x_{U},y_{U}]^{\rm
T}||_{2}}{|h_{H}-h_{U}|})})-{\sin}(\upsilon_{{\rm AoA}}))$, $\forall n$, in
the stage of designing the pilot matrix. Then, the pilot matrix to be
transmitted by the BS can be expressed as

${{\bm{U}_{0}^{\rm H}}={\bm{SDR}}({{{\bm{\bar{H}}}^{\rm H}}{{\bm{\Theta}}^{\rm
H}}\bm{A}(\bm{0})})^{-1}}$ (81)

###  -B Off-Grid Offsets Update

Given the estimated
$\hat{p}({\bm{x}})=\mathcal{CN}(0;{\bm{x}}^{post}_{B},v_{B}^{post}\bm{I})$ at
the $i$-th iteration, the progressively refined function in (49) can be
computed by

$\begin{array}[]{l}\hat{r}(\Delta\bm{\phi};\Delta{\bm{\phi}^{i}})\propto{{\rm{E}}_{\hat{p}(\bm{x}|\bm{y},\Delta\bm{\phi})}}[-\frac{1}{{\sigma_{e}^{2}}}||\bm{y}-\bm{F}(\Delta\bm{\phi})\bm{x}|{|^{2}}]\\\
\propto-\frac{1}{{\sigma_{e}^{2}}}(||\bm{y}-\bm{F}(\Delta\bm{\phi})\bm{x}_{B}^{post}|{|^{2}}+v_{B}^{post}{\rm{tr}}(\bm{F}(\Delta\bm{\phi})\bm{F}{(\Delta\bm{\phi})^{\rm
H}}))\end{array}$ (82)

The derivative of $\hat{r}(\Delta\bm{\phi};\Delta\bm{\phi}^{i})$ in (84) can
be calculated as
$\bm{\psi}^{(i)}_{\Delta\bm{\phi}}=[\psi^{(i)}(\Delta\phi_{1}),...,\psi^{(i)}(\Delta\phi_{N})]^{T}$,
where

$\begin{array}[]{l}{\psi^{(i)}}(\Delta{\phi_{n}})=2{\rm{Re}}({{\bm{a}}^{\prime}}{({{\hat{\omega}}_{n}}+\Delta{\phi_{n}})^{\rm
H}}{{\bm{U}}_{\bm{0}}}{{\bm{\bar{H}}}{\bm{\Theta}}}{({{\bm{U}}_{{0}}}{{\bm{\bar{H}}}{\bm{\Theta}}})^{\rm
H}}{\bm{a}}({{\hat{\omega}}_{n}}+\Delta{\phi_{n}}))\\\
\times\varphi_{1}^{(i)}+2{\rm{Re}}({{\bm{a}}^{\prime}}{({{\hat{\omega}}_{j}}+\Delta{\phi_{j}})^{\rm
H}}{{\bm{U}}_{{0}}}{{\bm{\bar{H}}}{\bm{\Theta}}}\bm{\varphi}_{2}^{(i)})\end{array}$
(83)

with
$\varphi_{1}^{(i)}=-\frac{1}{\sigma_{e}^{2}}(|{x}^{post}_{B,n}|^{2}+v_{B}^{post})$,
${\bm{\varphi}}_{2}^{(i)}=\frac{1}{\sigma_{e}^{2}}({x}^{post}_{B,n})*{\bm{y}}_{-n}$,
${\bm{y}}_{-n}={\bm{y}}-({\bm{U_{0}}}{\bm{\bar{H}}}{\bm{\Theta}})^{\rm
H}\sum_{j\neq
n}a(\widehat{\omega}_{j}+\Delta\phi_{j}){x}^{post}_{B,j},{\bm{a}}^{\prime}(\widehat{\omega}_{n}+\Delta\phi_{n})={\rm
d}{\bm{a}}(\widehat{\omega}_{n}+\Delta\phi_{n})/{\rm d}\Delta\phi_{n}$.

Then we can update the off-grid offsets as

$\displaystyle\Delta\bm{\phi}^{i+1}=\Delta\bm{\phi}^{i}+\alpha_{i}\bm{\psi}^{(i)}_{\Delta\bm{\phi}}$
(84)

###  -C Proof of Lemma 2

Define the matrix $\bm{\Lambda}=\bm{h}^{\rm H}\bm{\Theta}\bm{H}$. Then, we can
rewrite ${\rm SNR}_{u}$ in (3) as

${SNR}_{u}=P_{b}G\bm{v}^{\rm H}\bm{\Lambda}^{\rm
H}\bm{\Lambda}\bm{v}/\sigma_{0}^{2}$ (85)

According to (6) and (9), $\bm{\Lambda}^{\rm H}\bm{\Lambda}$ is a rank-1
martix with the following structure

$\bm{\Lambda}^{\rm H}\bm{\Lambda}=D\bm{a}_{\rm BS}(\upsilon_{{\rm
AoD}})\bm{a}_{\rm BS}(\upsilon_{{\rm AoD}})^{\rm H}$ (86)

where $D$ is a positive constant.

The non-zero eigenpair of $\bm{\Lambda}^{\rm H}\bm{\Lambda}$ can then be
expressed as

$\displaystyle(E_{val},E_{vec})=(D||\bm{a}_{\rm BS}(\upsilon_{{\rm
AoD}})||^{2}_{2},\frac{\bm{a}_{\rm BS}(\upsilon_{{\rm AoD}})}{||\bm{a}_{\rm
BS}(\upsilon_{{\rm AoD}})||_{2}})$ (87)

Owing to the unit magnitude of $\boldsymbol{v}$ and the Rayleigh-Ritz theorem,
$\boldsymbol{v}$ can be calculated by

$\displaystyle\bm{v}=\frac{\bm{a}_{\rm BS}(\upsilon_{{\rm
AoD}})}{||\bm{a}_{\rm BS}(\upsilon_{{\rm AoD}})||_{2}}$ (88)

###  -D Proof of Lemma 3

Define ${\mathcal{P}}$ as the feasible region of (70), and $\eta_{\rm
EE}=\min\limits_{k\in{\mathcal{K}}}\text{
}\displaystyle\frac{R_{u}(1-\epsilon_{u})+R_{k}(1-\epsilon_{k})}{{\bar{p}}_{b}+{\bar{p}}_{u}}=\min\limits_{k\in{\mathcal{K}}}\text{
}\frac{V_{k}(\bar{p}_{u})}{S(\bar{p}_{u})}$. Denote by $\eta_{\rm EE}^{\star}$
and $(\bar{p}_{u}^{\star}P_{U})\in{\mathcal{P}}$ the optimal energy efficiency
and the optimal UAV transmit power, respectively. Then, we have

$\displaystyle\eta_{\rm
EE}^{\star}=\displaystyle\frac{\min\limits_{k\in{\mathcal{K}}}V_{k}(\bar{p}_{u}^{\star})}{S(\bar{p}_{u}^{\star})}\geq\displaystyle\frac{\min\limits_{k\in{\mathcal{K}}}V_{k}(\bar{p}_{u})}{S(\bar{p}_{u})}$
(89)

Accordingly, the following equation holds

$\displaystyle\max\limits_{P_{u}}\text{
}{\min\limits_{k\in{\mathcal{K}}}\text{ }V_{k}(\bar{p}_{u})}-\eta_{\rm
EE}^{\star}{S(\bar{p}_{u})}=0$ (90)

On the contrary, suppose that $(\bar{p}_{u}^{\star}P_{U})$ is the optimal UAV
transmit power of the following optimization problem

$\displaystyle\max\limits_{P_{u}^{\prime}}\text{
}\min\limits_{k\in{\mathcal{K}}}\text{ }\displaystyle\text{
}V_{k}(P_{u}^{\prime}/{P_{U}})-\eta_{\rm EE}^{\star}S(P_{u}^{\prime}/{P_{U}})$
(91a) $\displaystyle{\rm s.t:}\text{ }{\rm Constraints\text{
}(\ref{eq_UAV_power_cons}),(\ref{eq:formualted_problem}c).}$ (91b)

and the following equation holds

$\displaystyle{\min\limits_{k\in{\mathcal{K}}}\text{
}V_{k}(\bar{p}_{u}^{\star})}-\eta_{\rm EE}^{\star}{S(\bar{p}_{u}^{\star})}=0$
(92)

As (70) and (91) have the same feasible region ${\mathcal{P}}$, for any UAV
transmit power $P_{u}^{\prime}\in{\mathcal{P}}$, the following inequality
holds

$\begin{array}[]{l}\mathop{\min}\limits_{k\in{\mathcal{K}}}\text{
}{\rm{}}{V_{k}}({P_{u}^{\prime}}/{P_{U}})-\eta_{{\rm{EE}}}^{\star}S({P_{u}^{\prime}}/{P_{U}})\leq\mathop{\min}\limits_{k\in{\mathcal{K}}}\text{
}{\rm{}}{V_{k}}(\bar{p}_{u}^{\star})-\eta_{{\rm{EE}}}^{\star}S(\bar{p}_{u}^{\star})=0\end{array}$
(93)

From (93), we obtain
$\eta_{{\rm{EE}}}^{\star}=\frac{{\mathop{\min}\limits_{k\in{\mathcal{K}}}{\rm{}}{V_{k}}(\bar{p}_{u}^{\star})}}{{S(\bar{p}_{u}^{\star})}}\geq\frac{{\mathop{\min}\limits_{k\in{\mathcal{K}}}{\rm{}}{V_{k}}(P_{u}^{\prime}/{P_{U}})}}{{S({P_{u}^{\prime}}/{P_{U}})}}$.
Then we can conclude that (70) and (91) have the same optimal UAV transmit
power. This completes the proof.

## References

  * [1] A. Brighente, J. Mohammadi, P. Baracca, S. Mandelli, and S. Tomasin, “Interference prediction for low-complexity link adaptation in beyond 5g ultra-reliable low-latency communications,” _IEEE Trans. Wirel. Commun._ , vol. 21, no. 10, pp. 8403–8415, 2022.
  * [2] G. K. Kurt, M. G. Khoshkholgh, S. Alfattani, A. Ibrahim, T. S. J. Darwish, M. S. Alam, H. Yanikomeroglu, and A. Yongaçoglu, “A vision and framework for the high altitude platform station (HAPS) networks of the future,” _IEEE Commun. Surv. Tutorials_ , vol. 23, no. 2, pp. 729–779, 2021.
  * [3] 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 Trans. Ind. Informatics_ , vol. 16, no. 9, pp. 5866–5873, 2020.
  * [4] N. S. Perovic, L. Tran, M. D. Renzo, and M. F. Flanagan, “Achievable rate optimization for MIMO systems with reconfigurable intelligent surfaces,” _IEEE Trans. Wirel. Commun._ , vol. 20, no. 6, pp. 3865–3882, 2021.
  * [5] A. Ranjha, G. Kaddoum, and K. Dev, “Facilitating URLLC in UAV-assisted relay systems with multiple-mobile robots for 6G networks: A prospective of agriculture 4.0,” _IEEE Trans. Ind. Informatics_ , vol. 18, no. 7, pp. 4954–4965, 2022.
  * [6] C. M. W. Basnayaka, D. N. K. Jayakody, and Z. Chang, “Age-of-information-based URLLC-enabled UAV wireless communications system,” _IEEE Internet Things J._ , vol. 9, no. 12, pp. 10 212–10 223, 2022.
  * [7] G. Bartoli and D. Marabissi, “CQI prediction through recurrent neural network for UAV control information exchange under URLLC regime,” _IEEE Trans. Veh. Technol._ , vol. 71, no. 5, pp. 5101–5110, 2022.
  * [8] K. Wang, C. Pan, H. Ren, W. Xu, L. Zhang, and A. Nallanathan, “Packet error probability and effective throughput for ultra-reliable and low-latency UAV communications,” _IEEE Trans. Commun._ , vol. 69, no. 1, pp. 73–84, 2021\.
  * [9] C. She, C. Liu, T. Q. S. Quek, C. Yang, and Y. Li, “Uav-assisted uplink transmission for ultra-reliable and low-latency communications,” in _ICC Workshops_. IEEE, 2018, pp. 1–6.
  * [10] P. Yang, X. Xi, K. Guo, T. Q. S. Quek, J. Chen, and X. Cao, “Proactive UAV network slicing for URLLC and mobile broadband service multiplexing,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 10, pp. 3225–3244, 2021.
  * [11] A. Ranjha and G. Kaddoum, “Quasi-optimization of uplink power for enabling green URLLC in mobile UAV-assisted IoT networks: A perturbation-based approach,” _IEEE Internet Things J._ , vol. 8, no. 3, pp. 1674–1686, 2021\.
  * [12] K. Chen, Y. Wang, J. Zhao, X. Wang, and Z. Fei, “URLLC-oriented joint power control and resource allocation in UAV-assisted networks,” _IEEE Internet Things J._ , vol. 8, no. 12, pp. 10 103–10 116, 2021.
  * [13] X. Xi, X. Cao, P. Yang, J. Chen, T. Q. S. Quek, and D. Wu, “Network resource allocation for eMBB payload and URLLC control information communication multiplexing in a multi-uav relay network,” _IEEE Trans. Commun._ , vol. 69, no. 3, pp. 1802–1817, 2021.
  * [14] P. Yang, X. Xi, T. Q. S. Quek, J. Chen, and X. Cao, “Power control for a URLLC-enabled UAV system incorporated with DNN-based channel estimation,” _IEEE Wirel. Commun. Lett._ , vol. 10, no. 5, pp. 1018–1022, 2021.
  * [15] Y. Li, C. Yin, T. Do-Duy, A. Masaracchia, and T. Q. Duong, “Aerial reconfigurable intelligent surface-enabled URLLC UAV systems,” _IEEE Access_ , vol. 9, pp. 140 248–140 257, 2021.
  * [16] A. Ranjha and G. Kaddoum, “URLLC facilitated by mobile UAV relay and RIS: A joint design of passive beamforming, blocklength, and UAV positioning,” _IEEE Internet Things J._ , vol. 8, no. 6, pp. 4618–4627, 2021.
  * [17] X. Zhang, J. Wang, and H. V. Poor, “Statistical QoS-driven beamforming and trajectory optimizations in UAV/IRS-based 6G wireless networks in the non-asymptotic regime,” in _ISIT_. IEEE, 2022, pp. 3333–3338.
  * [18] M. Samir, M. K. Elhattab, C. Assi, S. Sharafeddine, and A. Ghrayeb, “Optimizing age of information through aerial reconfigurable intelligent surfaces: A deep reinforcement learning approach,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 4, pp. 3978–3983, 2021.
  * [19] H. Jeon, S. Park, J. Park, K. Huang, and C. Chae, “An energy-efficient aerial backhaul system with reconfigurable intelligent surface,” _IEEE Trans. Wirel. Commun._ , vol. 21, no. 8, pp. 6478–6494, 2022.
  * [20] N. Gao, S. Jin, X. Li, and M. Matthaiou, “Aerial RIS-assisted high altitude platform communications,” _IEEE Wirel. Commun. Lett._ , vol. 10, no. 10, pp. 2096–2100, 2021.
  * [21] A. Khalili, E. M. Monfared, S. Zargari, M. R. Javan, N. M. Yamchi, and E. A. Jorswieck, “Resource management for transmit power minimization in UAV-assisted RIS HetNets supported by dual connectivity,” _IEEE Trans. Wirel. Commun._ , vol. 21, no. 3, pp. 1806–1822, 2022.
  * [22] X. Hu, R. Zhang, and C. Zhong, “Semi-passive elements assisted channel estimation for intelligent reflecting surface-aided communications,” _IEEE Trans. Wirel. Commun._ , vol. 21, no. 2, pp. 1132–1142, 2022.
  * [23] S. Li, B. Duo, M. D. Renzo, M. Tao, and X. Yuan, “Robust secure UAV communications with the aid of reconfigurable intelligent surfaces,” _IEEE Trans. Wirel. Commun._ , vol. 20, no. 10, pp. 6402–6417, 2021.
  * [24] D. Grace and M. Mohorcic, _Broadband communications via high altitude platforms_. John Wiley & Sons, 2011.
  * [25] D. Tse and P. Viswanath, _Fundamentals of Wireless Communication_. Cambridge University Press, 2005.
  * [26] A. Al-Hourani, K. Sithamparanathan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” _IEEE Wirel. Commun. Lett._ , vol. 3, no. 6, pp. 569–572, 2014.
  * [27] Y. Polyanskiy, H. V. Poor, and S. Verdú, “Channel coding rate in the finite blocklength regime,” _IEEE Trans. Inf. Theory_ , vol. 56, no. 5, pp. 2307–2359, 2010.
  * [28] W. Yang, G. Durisi, T. Koch, and Y. Polyanskiy, “Quasi-static multiple-antenna fading channels at finite blocklength,” _IEEE Transactions on Information Theory_ , vol. 60, no. 7, pp. 4232–4265, 2014.
  * [29] L. Yuan, N. Yang, F. Fang, and Z. Ding, “Performance analysis of UAV-assisted short-packet cooperative communications,” _IEEE Trans. Veh. Technol._ , vol. 71, no. 4, pp. 4471–4476, 2022.
  * [30] C. Pan, H. Ren, Y. Deng, M. Elkashlan, and A. Nallanathan, “Joint blocklength and location optimization for URLLC-enabled UAV relay systems,” _IEEE Commun. Lett._ , vol. 23, no. 3, pp. 498–501, 2019.
  * [31] J. Cheng, C. Shen, Z. Chen, and N. Pappas, “Robust beamforming design for IRS-aided URLLC in D2D networks,” _IEEE Trans. Commun._ , vol. 70, no. 9, pp. 6035–6049, 2022.
  * [32] R. Hashemi, S. Ali, N. H. Mahmood, and M. Latva-aho, “Joint sum rate and blocklength optimization in RIS-aided short packet URLLC systems,” _IEEE Commun. Lett._ , vol. 26, no. 8, pp. 1838–1842, 2022.
  * [33] W. Chang, Z. Meng, K. Liu, and L. Wang, “Energy-efficient sleep strategy for the UBS-assisted small-cell network,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 5, pp. 5178–5183, 2021.
  * [34] Z. Wang, R. Liu, Q. Liu, J. S. Thompson, and M. Kadoch, “Energy-efficient data collection and device positioning in UAV-assisted IoT,” _IEEE Internet Things J._ , vol. 7, no. 2, pp. 1122–1139, 2020.
  * [35] J. Chakareski, S. Naqvi, N. Mastronarde, J. Xu, F. Afghah, and A. Razi, “An energy efficient framework for UAV-assisted millimeter wave 5G heterogeneous cellular networks,” _IEEE Trans. Green Commun. Netw._ , vol. 3, no. 1, pp. 37–44, 2019.
  * [36] K. Ardah, S. Gherekhloo, A. L. F. de Almeida, and M. Haardt, “TRICE: A channel estimation framework for RIS-aided millimeter-wave MIMO systems,” _IEEE Signal Process. Lett._ , vol. 28, pp. 513–517, 2021.
  * [37] P. Wang, J. Fang, H. Duan, and H. Li, “Compressed channel estimation for intelligent reflecting surface-assisted millimeter wave systems,” _IEEE Signal Process. Lett._ , vol. 27, pp. 905–909, 2020.
  * [38] A. Taha, M. Alrabeiah, and A. Alkhateeb, “Enabling large intelligent surfaces with compressive sensing and deep learning,” _IEEE Access_ , vol. 9, pp. 44 304–44 321, 2021.
  * [39] P. An, X. Cao, P. Yang, K. Guo, Y. Gao, and T. Q. S. Quek, “Energy-efficient URLLC service provision via a near-space information network,” Tech. Rep. 2304.04163v2, 2023, https://arxiv.org/pdf/2304.04163.pdf.
  * [40] S. M. Kay, _Fundamentals of statistical signal processing: estimation theory_. Prentice-Hall, Inc., 1993.
  * [41] J. Ma and L. Ping, “Orthogonal AMP,” _IEEE Access_ , vol. 5, pp. 2020–2033, 2017.
  * [42] J. Dai, A. Liu, and V. K. N. Lau, “FDD massive MIMO channel estimation with arbitrary 2D-array geometry,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2584–2599, 2018.
  * [43] C. J. Wu, “On the convergence properties of the EM algorithm,” _The Annals of statistics_ , pp. 95–103, 1983.
  * [44] S. Dhok, P. Raut, P. K. Sharma, K. Singh, and C. Li, “Non-linear energy harvesting in RIS-assisted URLLC networks for industry automation,” _IEEE Trans. Commun._ , vol. 69, no. 11, pp. 7761–7774, 2021.
  * [45] Y. Gu, H. Chen, Y. Li, and B. Vucetic, “Ultra-reliable short-packet communications: Half-duplex or full-duplex relaying?” _IEEE Wirel. Commun. Lett._ , vol. 7, no. 3, pp. 348–351, 2018.
  * [46] Y. Ye, _Interior point algorithms - theory and analysis_ , ser. Wiley-Interscience series in discrete mathematics and optimization. Wiley, 1998.
  * [47] G. Z. Karabulut and A. Yongacoglu, “Sparse channel estimation using orthogonal matching pursuit algorithm,” in _IEEE 60th Vehicular Technology Conference, 2004. VTC2004-Fall. 2004_ , vol. 6. IEEE, 2004, pp. 3880–3884.
  * [48] W. Dai and O. Milenkovic, “Subspace pursuit for compressive sensing signal reconstruction,” _IEEE Trans. Inf. Theory_ , vol. 55, no. 5, pp. 2230–2249, 2009.
  * [49] A. Alsharoa and M. Alouini, “Improvement of the global connectivity using integrated satellite-airborne-terrestrial networks with resource optimization,” _IEEE Trans. Wirel. Commun._ , vol. 19, no. 8, pp. 5088–5100, 2020.
  * [50] J. Ziniel and P. Schniter, “Dynamic compressive sensing of time-varying signals via approximate message passing,” _IEEE Trans. Signal Process._ , vol. 61, no. 21, pp. 5270–5284, 2013.