Unmanned Aerial Vehicle-Enabled Aerial Radio Environment Map Construction: A Multi-Stage Approach to Data Sampling and Path Planning

Abstract

An aerial Radio Environment Map (REM) characterizes the spatial distribution of Received Signal Strength (RSS) across a geographic space of interest, which is crucial for optimizing wireless communication in the air. Aerial REM construction can rely on Unmanned Aerial Vehicles (UAVs) to autonomously select interesting positions for sampling RSS data, enhancing the quality of construction. However, due to the lack of prior information about the environment, it is challenging for UAVs to determine suitable sampling positions online. Additionally, achieving efficient exploration of the target area through collaboration among multiple UAVs is difficult. To address this issue, this paper proposes a multi-stage approach to data sampling and path planning with multiple UAVs. Specifically, the UAVs’ data sampling task over the target area is divided into multiple stages. By selecting an appropriate stage position, we use the RSS values at that position to determine whether additional data need to be sampled in a specific local area. At each stage, the area is divided into Voronoi diagrams based on the current position of each UAV, assigning each UAV its own region to explore. In our sampling strategy, the probability distribution for sampling is obtained by estimating the RSS and uncertainty of unsampled positions and then taking the weighted sum of these two values. To obtain the shortest flight path for selected sampling positions, we employ a network structure based on self-attention as the policy network, which is trained through the actor–critic framework to obtain an improvement heuristic strategy, replacing traditional manually designed strategies. Experimental results across three different scenarios indicate that the approach improves the quality of aerial REM construction while efficiently planning the shortest paths for UAVs between sampling positions.

1. Introduction

A Radio Environment Map (REM) characterizes the spatial distribution of the Received Signal Strength (RSS) across a geographic space of interest [1]. It has been widely used for communication network optimization, such as the deployment of base stations [2], reducing communication interference [3], and user localization [4]. With the continuous expansion of human activities towards the sky and the usage of Unmanned Aerial Vehicles (UAVs) for various applications [5], constructing an aerial REM helps improve the communication service quality in the sky. Some research works have been conducted on how to use an aerial REM to optimize communication between aerial devices and ground-based stations [5,6,7].

A typical REM is constructed by ground-based sensors with measured RSS data, which are not suitable for aerial REM construction. Fortunately, UAVs equipped with Spectrum Monitoring Devices (SMDs) can be used to sample the RSS data in air space, enabling the construction of a high-quality aerial REM. UAVs can dynamically plan flight paths, selecting the most informative measurement positions for data sampling. In the task of UAV-assisted RSS data sampling, most methods assume the partial prior information of the environment is known. For instance, in the work [8], it is assumed that the locations of the radiation sources are known, and then a data-driven deep learning algorithm is utilized to measure the uncertainty at positions within the target area, indicating the informative measurement at each position for UAV selection. The authors in [9] predetermine some areas as regions of interest, assigning higher informative value to the positions within these regions. When there is a lack of prior information, most work involves sampling by randomly selecting positions [10] or following a fixed sampling trajectory [11].

However, the radiation source is unknown in most applications [12], and UAVs may struggle to determine the most informative measurement positions in an online way. Additionally, some studies do not propose a path planning method to efficiently sample data after optimizing the sampling positions. Moreover, the above works were completed by a single UAV for data sampling, without considering the scenario of multi-UAV collaboration.

To address the above challenges, in this paper, we propose a multi-stage approach to data sampling and path planning with multiple UAVs. For effective collaboration among UAVs, we employ a Voronoi partitioning method. For the challenge of unknown prior environmental information, we propose a granular sampling strategy for selecting sampling positions. This selection is based on using RSS data sampled during flight to estimate the RSS value and uncertainty at other unsampled positions. For the problem of planning the shortest flight path to the selected sampling positions, we use a self-attention-based network structure as a policy network to learn a heuristic strategy for path planning. An earlier version of this paper was accepted as a conference paper in ICCT 2024 [13]. This version is a major extension with focuses on the generalization of the algorithm across different datasets.

The contributions of this work are as follows:

• A Voronoi partitioning method is proposed to efficiently coordinate multiple UAVs for data sampling in the target area. Each UAV uses the Voronoi partitioning to focus on different local regions to avoid overlapping exploration areas.
• A granular sampling strategy is proposed. It utilizes RSS data sampled during flight to estimate the RSS and uncertainty at unsampled locations. The weighted sum of both is used to obtain the probability distribution for sampling, ensuring that sampling positions are distributed as widely as possible within the areas of interest.
• A network architecture based on a self-attention mechanism is proposed as a policy network, which is trained offline using the actor–critic framework to obtain an improved heuristic strategy, replacing traditional manually designed strategies.
• Extensive simulations show that the proposed approach facilitates cooperation among multiple UAVs, enhances the construction accuracy of aerial REMs, and enables the planning of shorter flight paths for UAVs.

2. Related Works

2.1. REM Construction

The methods for constructing REMs are primarily divided into two categories: model-driven and data-driven approaches. Model-driven methods are based on physics simulations, calculating wireless channels using Maxwell’s equations. By inputting the 3D model of the environment along with the positions of transmitters and receivers, software can generate the REM of the input environment [14]. However, model-driven methods have high computational complexity and require a significant amount of prior environmental information. On the contrary, data-driven methods do not require prior environmental information; they construct the REM directly from sampling data. In recent years, geological interpolation algorithms, such as the inverse distance weighting (IDW) algorithm [15], Gaussian Process Regression (GPR) [16], and the kriging algorithm (KGA) [17], have been widely used in the construction of REMs. With the advancement of artificial intelligence, there have been some recent works utilizing data-driven methods based on machine learning to construct REMs, such as the generative adversarial network (GAN) [18], the graph neural network (GNN) [19], the autoencoder [20], and the convolutional neural network (CNN) [21]. However, machine learning-based methods rely heavily on the quality of the dataset, and their robustness can be compromised when the environment changes.

In addition to various REM construction methods, some studies have focused on optimizing the distribution of sampling positions to improve the accuracy of REM construction. Some methods optimize the sampling positions by utilizing partial prior information of the environment [8,9,22,23,24]. For example, by assuming the position and power of the radiation source are known, a Bayesian estimation-based deep learning algorithm is used in [8] to estimate the uncertainty of each sampling position. The authors in [9] predefine certain areas as regions of interest and distributed the sampling positions within these regions as much as possible. In [22,23], the authors assume that some RSS data are collected in advance in the target area. In the work [24], the REM constructed from the previous time step is used as prior information. However, it is difficult to obtain the prior information in practice. When in an unknown environment, some work involves sampling through random sampling or by setting a fixed flight trajectory [1,10,11,25], while other work involves selecting sampling positions based on certain assumed conditions [26,27,28,29]. In [26], the authors assume an exponential decay model with respect to distance to represent the characteristics of the propagation channel, thereby determining suitable sampling positions. The authors in [27,28,29] model the distribution of the data to be sampled within the target area as a Gaussian process. However, due to interference between multiple radiation sources and the complexity of the propagation environment, this assumption is often difficult to establish.

2.2. UAV Trajectory Design

When using the sampling strategy proposed in this paper to determine the sampling points for UAVs, the trajectory planning problem of the UAV can be transformed into a path planning problem, similar to the classic Traveling Salesman Problem (TSP). Due to the NP-hard nature of the path planning problem [30], solving it presents a significant challenge. Classic approaches can be categorized into approximation methods and heuristic methods. Although approximation methods can quickly find solutions, the obtained solutions may deviate significantly from the optimal solution, especially for certain problems where the approximation quality can be poor [31]. Heuristic methods find local optimal solutions within a reasonable computation time by selecting a heuristic strategy [32]. However, the effectiveness of traditional heuristic algorithms depends on the chosen heuristic strategy and initial solution. Recently, there has been a growing trend toward applying deep learning to automatically discover heuristic algorithms for solving routing problems [33]. Through the learning capabilities of deep networks, Deep Reinforcement Learning (DRL) algorithms can be generated to solve problems using heuristic strategies instead of those designed by humans, such as the recurrent neural network (RNN) [34] and the attention mechanism [35].

2.3. Motivation

In summary, the use of UAVs to assist in REM construction presents various challenges. In addition to the selection of sampling positions and the design of UAV trajectory, issues such as multi-UAV communication synchronization, environmental interference, and adaptability to complex terrain must also be addressed. However, many existing works have explored relevant solutions. Regarding multi-UAV communication synchronization, Jin et al. [36] proposed a scheme for assisting time synchronization in multi-UAV networks by utilizing frequency offset information. In terms of environmental interference, Li et al. [37] developed an alternating optimization approach to ensure reliable communication for UAVs in complex environments with low signal-to-noise ratios. The authors in [38] also investigated the impact of adverse weather conditions on UAV communication channels. Concerning adaptability to complex terrain, Zhang et al. [39] proposed an approach that combines reinforcement learning and artificial potential field methods, enabling UAVs to achieve safe flight in dynamic and unknown environments. Unlike these studies, the focus of our work is on addressing the challenges of coordination among multiple UAVs, selecting appropriate sampling positions, and designing suitable trajectories for the UAVs. In existing related works, the sampling positions of UAVs are either random or fixed, or optimized based on prior information or assumptions, without considering the possibility of dynamically optimizing the sampling positions using RSS data sampled by the UAVs during flight. Moreover, when the sampling locations are determined, many studies rely on manually designed heuristic strategies to plan flight trajectories. Not only do these approaches result in long computation times, but the solution quality also depends on the initial solution and the designed strategy. To the best of our knowledge, none of the studies have considered scenarios involving multiple UAVs assisting in REM construction while simultaneously addressing both the sampling position optimization and the path planning problem. To fill these gaps, this paper proposes a multi-stage approach to data sampling and path planning. It achieves multi-UAV coordination through Voronoi partitioning, optimizes sampling positions based on uncertainty, and designs a policy network based on self-attention. Using the actor–critic training method, it learns a heuristic strategy offline to quickly plan the shortest path between UAV sampling positions.

3. System Model

We consider the urban scenario with multiple-UAV-enabled REM construction, as depicted in Figure 1. There are a total of U UAVs and S radio sources, located within a square region with a side length of Q. Let $\mathcal{U}$ = {1, ..., U} and $\mathcal{S}$ = {1, ..., S} denote the sets of UAVs and radiation sources, respectively. We divide the target region into small square grids with side length l, resulting in a total of $D\times D$ grids within the region, where $D=Q/l$. In order to construct the aerial REM, UAVs equipped with SMD collect the RSS data at the centers of small square grids.

Table 1. List of notation.

Symbol	Definition	Symbol	Definition
$\mathcal{U},\mathcal{S},{\mathcal{G}}_D$	Sets of UAVs, radio sources, and grids, respectively.	$Q,l$	Side length of the region and small grids, respectively.
${\mathbf{G}}_d(i,j)$	Position of the (i, j)-th grid.	$U,H$	Number and altitude of UAVs, respectively.
${\Omega }^u$	Set of collecting positions for UAV u.	$P_{ij}$	RSS value at the ($i_s$, $j_s$)-th position in the grid.
${\mathbf{a}}^u$	Vector of the collecting selection of UAV u.	$\mathbf{X}$	Set of measurement data points by all UAVs.

lists the notations in this paper.

3.1. REM Construction Model

The grid for discrete positions of the region area can be represented as ${\mathcal{G}}_D=\{{\mathbf{G}}_d(i,j)\in {\mathbb{R}}^{2\times 1}:i,j\in \mathcal{D}\}$, where $\mathcal{D}=\{1,\cdots ,D\}$. Here, ${\mathbf{G}}_d(i,j)$ represents the position of the (i, j)-th grid and is expressed by the following equation:

$${\mathbf{G}}_d(i,j)={[i,j]}^{\mathrm{T}}l,i,j\in \mathcal{D}.$$

(1)

As for UAV u, it moves horizontally at a fixed altitude H. Denote ${\Omega }^u$ as the set of sampling positions for UAV u, and the path length of UAV u is denoted as ${\sum }_{k=1}^{\left|{\Omega }^u\right|}∥{\Omega }_k^u-{\Omega }_{k-1}^u∥$, where ${\Omega }_k^u$ represents UAV u sampling at the k-th sampling location. Furthermore, denote ${\mathbf{a}}^u=\{a_{ij}^u:i,j\in \mathcal{D}\}$ as the vector of the sampling selection of UAV u, where $a_{ij}^u=1$ denotes that UAV u samples the RSS at ${\mathbf{G}}_d(i,j)$, which means that ${\mathbf{G}}_d(i,j)\in {\Omega }^u$ and $a_{ij}^u=0$ otherwise. The RSS value sampled by the UAV at ${\mathbf{G}}_d(i,j)$ is defined as $P_{ij}$.

The principle of the method for constructing an REM is to find an optimal function $f(·)$, using the set $\mathbf{X}$ of measurement data points, where $\mathbf{X}=\{{[il,jl,P_{ij}]}^{\mathrm{T}}\mid a_{ij}^u=1,\exists u\in \mathcal{U}\}$, to construct the REM $\tilde{p}\in {\mathbb{R}}^{D\times D}$. The overall process of REM construction can be expressed as min ${∥\tilde{p}-p∥}_F$, where $p\in {\mathbb{R}}^{D\times D}$ is the real REM and ${∥·∥}_F$ means the Frobenius norm for matrices.

3.2. Problem Formulation

Based on the established system model, multiple-UAV-enabled REM construction mainly focuses on two aspects, i.e., how to find appropriate sampling positions for UAVs to improve the quality of aerial REM construction and how to find a shortest path when determining the positions of sampling points for UAVs. These goals are achieved by optimizing the positions and sequence of UAV sampling points. The problem can be expressed as:

$$\underset{{\mathbf{a}}^u,{\Omega }^u}{min}\sigma \sum _{u=1}^U\sum _{k=1}^{\left|{\Omega }^u\right|}∥{\Omega }_k^u-{\Omega }_{k-1}^u∥+\tau ∥\tilde{p}-p∥F,$$

(2)

$$\mathrm{s}.\mathrm{t}.\tilde{p}=f\left(\mathbf{X}\right),$$

(3)

$$\mathrm{s}.\mathrm{t}.\sum {u=1}^U\sum _{i=1}^D\sum _{j=1}^D{a}_{ij}^u\le \epsilon D^2,$$

(4)

In Equation (2), the first term aims to minimize the UAV flight path while optimizing the sampling positions. The second item focuses on minimizing the discrepancy between the constructed REM and the real REM by optimizing the selection vector ${\mathbf{a}}^u$. Constraint (3) stipulates that the construction of the REM must utilize data collected from all UAVs, as processed by the function $f(·)$. Constraint (4) states that the total volume of data sampled by all UAVs must not exceed the sampling rate $\epsilon$, reflecting the practical limitation that not every position in the target area can be sampled.

In fact, Equation (2) represents a multi-objective optimization challenge. It is particularly difficult to determine optimal sampling positions without knowing the locations of radiation sources. Moreover, even after identifying these sampling positions, planning the shortest route between them is an NP-hard problem. To tackle these issues, we propose a multi-stage approach to data sampling and path planning that utilizes only the RSS measurements and the locations of UAVs to identify suitable sampling positions. Once these positions are established, we employ a DRL approach to design the most efficient flight path.

4. Multi-Stage Approach to Data Sampling and Path Planning

In the task of aerial REM construction in an unknown environment, data-driven methods often emphasize the spatial characteristics of the data [40], which usually implies that the RSS values exhibit a certain degree of spatial continuity on the map. However, various factors during signal propagation can disrupt the continuity, such as interference between multiple radiation sources and obstruction by obstacles [41].

To improve the construction quality of the REM, we can adopt a granular sampling strategy in areas where RSS values change dramatically, meaning sampling more data in these areas. This can be effectively achieved by dividing the data sampling task with UAVs for the target area into multiple stages. Multi-stage refers to the process of selecting several appropriate stage positions for each UAV. In this way, the entire target area is divided into multiple local areas. Each UAV can determine whether a granular sampling strategy is needed for the local area corresponding to the current stage based on the change in RSS values between the two consecutive stage positions. The steps in each stage are as follows:

• At the beginning of each stage s, based on the current positions $O_s$ of each UAV, Voronoi partitioning is applied to assign a region of responsibility to each UAV, restricting it to select the starting position $O_{s+1}$ for the next stage within its partitioned region.
• When UAV u flies from $O_s$ to $O_{s+1}$, it will decide whether to adopt a granular sampling strategy based on the magnitude of the change in RSS values between its current position $O_s$ and its previous position $O_{s-1}$. The granular sampling strategy uses the already-sampled RSS data to construct the RSS values and uncertainties at unsampled positions. The two values are weighted and summed to obtain a probability distribution for sampling. Then, UAV u selects appropriate sampling positions based on the probability distribution.
• After determining the next stage position $O_{s+1}$ and sampling positions for UAV u at stage s, a DRL method is proposed for planning the shortest path at stage s. We first randomly connect these positions into a path and then train our strategy using the actor–critic algorithm. The trained strategy is used to perform 2-opt optimization on the randomly generated path, ultimately planning the shortest path for the UAV u to fly from $O_s$ to $O_{s+1}$ and complete the sampling task at stage s.

4.1. Stage Position Selection Based on Voronoi Partitioning

In selecting stage positions, we want the UAV to choose an appropriate position to determine whether the local area in that stage requires additional sampling. However, if multiple UAVs are present, they might move in the same direction, exploring overlapping areas, which can lead to inefficiency and longer task completion times. To avoid this, we use a Voronoi partitioning [42] approach. For any given UAV u, a Voronoi partition is generated based on its own position, and it only considers boundary positions within its own Voronoi partition for the next move. This will effectively prevent the UAVs from exploring the same areas repeatedly.

For UAV u, each stage s focuses on a circular area centered at $O_s$ with radius r, where r depends on the length Q of the mission region. UAV u will select the starting position $O_{s+1}$ for the next stage along the edge of the circle. We partition the area into multiple sections using Voronoi partitioning at the beginning of each stage s. Each UAV is restricted to selecting its stage position $O_{s+1}$ from within its own partitioned area. Additionally, to prevent the UAVs from revisiting previously explored regions, we designate these areas as explored in subsequent stages. The partitioning of the UAVs is depicted in Figure 2.

By discretizing the reachable positions of the entire area for the UAV, UAV u can obtain a coordinate set (${\mathrm{X}}_{ti}^u$,${\mathrm{Y}}_{ti}^u$) for the selectable stage positions at each stage s. To select the most suitable stage position from the set, we define the following utility function $F_t^u$:

$$F_t^u=\alpha d_{int}+(1-\alpha )d_n,$$

(5)

where $d_{int}$ denotes the distance between the candidate target position and the initial position of UAV u (i.e., the position of the first stage) and $d_n$ denotes the distance between the candidate target position and the nearest neighboring UAV position. Note that when $\alpha$ = 1 , the UAV performs a depth-first exploration, allowing it to explore deeper into the area. When $\alpha$ = 0, the UAV aims to maximize its distance from other UAVs, which helps prevent repeated exploration of the same area. Regarding the task environment, we set an appropriate value for $\alpha$, and the UAV iteratively searches for the position with the minimum $F_t^u$ as the next stage position $O_{s+1}$. When the entire area is marked as explored, it indicates that the task is completed, and UAVs conclude the cooperative data sampling mission. The stage position selection based on Voronoi partitioning is written in Algorithm 1.

Algorithm 1 Stage Position Selection Based on Voronoi Partitioning

1: Input: Current positions $O_s$ of each UAV. 2: Output: Next stage positions $O_{s+1}$ of each UAV. 3: while the entire area is not marked as explored do 4:  Partition via Voronoi diagrams based on the current positions $O_s$ of each UAV u. 5:  for each UAV $u\in \{1,\cdots ,U\}$ do 6:   Under the constraint of radius r, obtain the set of selectable next stage positions (${\mathrm{X}}_{ti}^u$,${\mathrm{Y}}_{ti}^u$) for UAV u within its own partition and the unexplored area; 7:   Calculate $F_t^u$ of each coordinate in the set (${\mathrm{X}}_{ti}^u$,${\mathrm{Y}}_{ti}^u$) through (5); 8:   Select the position with the minimum $F_t^u$ as the position $O_{s+1}$ for UAV u. 9:  end for 10: end while

4.2. The Granular Sampling Strategy Based on Uncertainty

The decision of whether to adopt a granular sampling strategy in stage s is based on the magnitude of the change in RSS values between its current position $O_s$ and its previous position $O_{s-1}$, which can be defined as:

$$\left\{\begin{array}{cc}\mathrm{Directly}\mathrm{fly}\mathrm{to}{O}_{s+1},\hfill & \mathrm{if}|R_s^u-R_{s-1}^u|<R_{th}\hfill \\ \begin{array}{c}\mathrm{Adopt}\mathrm{a}\mathrm{granular}\hfill \\ \mathrm{collecting}\mathrm{strategy},\hfill \end{array}\hfill & \mathrm{if}|R_s^u-R_{s-1}^u|\ge R_{th}\hfill \end{array}\right.$$

(6)

where $R_s^u$ and $R_{s-1}^u$ are the RSS values of UAV u at positions $O_s$ and $O_{s-1}$ and $R_{th}$ is the set threshold. When the magnitude of the change in RSS values is greater than $R_{th}$, UAV u will sample $n_{set}$ points in the current local area before moving to the next stage position $O_{s+1}$, where $n_{set}$ depends on the sampling rate $\epsilon$. The local area is a square region with the stage position $O_s$ as the center and side length $2r$.

In our multi-UAV system, we assume that each UAV is capable of broadcasting its stored point information and current position, as well as receiving broadcasts from other UAVs. Based on the above assumptions, there may be $n_s^u$ already-sampled points within the local area of UAV u at stage s. To choose suitable sampling positions, we will use the sampled points to construct the REM of the local area by the KGA, not only constructing the REM but also representing the uncertainty of each position through variance values. The KGA can be expressed as:

$$\widehat{P_0}=\sum _{n=1}^{n_s^u}{\lambda }_n{P}_n,$$

(7)

where $\widehat{P_0}$ is the estimated RSS value at the unsampled position $(x_0,y_0)$ and ${\lambda }_n$ is the weight coefficient of the sampled position $(x_n,y_n)$. Due to the probability of establishing LoS between the UAV and ground radiation sources being greater in high altitude, we assume that for the RSS at any position in space, there is the same expected value c and variance ${\sigma }^2$. Then, we have:

$$P_n=c+R\left(n\right),$$

(8)

where $R\left(n\right)$ is the random deviation at point $(x_n,y_n)$, satisfying $\mathrm{Var}\left[R\left(n\right)\right]={\sigma }^2$. It is important to note that in real-world environments, LoS conditions may not hold in dense urban areas with high interference. When LoS conditions fail, additional losses due to shadowing and gains from small-scale fading and unmodeled effects must be considered, alongside $R\left(n\right)$ in Equation (8). In the KGA implementation, these various losses are modeled as a noise. As a result, KGA’s performance in dense urban environments with high interference will be worse compared to environments with lower interference. To overcome this limitation, it is necessary to combine model-based methods with more precise interference modeling of the environment or to rely on ray tracing methods to obtain more accurate numerical results.

Based on the assumptions stated above, the optimal weight coefficient should satisfy:

$$\left[\begin{array}{cccc}{r}_{11}& \cdots & r_{1n_s^u}& 1\\ ⋮& \ddots & ⋮& ⋮\\ r_{n_s^{u}1}& \cdots & r_{n_s^u{n}_s^u}& 1\\ 1& \cdots & 1& 0\end{array}\right]\left[\begin{array}{c}{\lambda }_1\\ ⋮\\ {\lambda }_{n_s^u}\\ -\varphi \end{array}\right]=\left[\begin{array}{c}{r}_{10}\\ ⋮\\ r_{n_s^{u}0}\\ 1\end{array}\right],$$

(9)

where $\varphi$ is the Lagrange multiplier and $r_{ij}$ is the semivariance function value between points $(x_i,y_i)$ and $(x_j,y_j)$, defined as:

$$r_{ij}={\displaystyle \frac{1}{2}}\mathbb{E}\left[{(P_i-P_j)}^2\right].$$

(10)

Since the RSS value at $(x_0,y_0)$ is unknown, the key to computing the weight ${\lambda }_n$ is to solve for ${[r_{10},\cdots ,r_{n_s^{u}0}]}^T$. In general, the RSS values between two positions in space are a function of the distance d between the two positions, which requires finding an optimal fitting curve to fit the relationship between d and r. Similar to [17], we assume that the shadow loss at each position is spatially correlated. When the shadow values at two positions are $W_i$ and $W_j$, respectively, the correlation coefficient ${\rho }_{ij}$ between them can be expressed as:

$${\rho }_{ij}={\displaystyle \frac{\mathbb{E}\left[W_i{W}_j\right]}{{\sigma }^2}}=\exp (-{\displaystyle \frac{d_{ij}}{\delta }}),$$

(11)

where $d_{ij}$ is the distance between the two positions, ${\sigma }^2$ is a constant, and $\delta$ is the distance at which the correlation decays to 1/2. Because of the assumption of correlated shadowing expressed in (11), an exponential semivariogram model is used for fitting:

$$r_{ij}=c_0+c_1(1-\exp (-{\displaystyle \frac{d_{ij}}{a}})),$$

(12)

where $c_0$, $c_1$, and a are the parameters that need to be fitted.

When we derive the weights ${[{\lambda }_1,\cdots ,{\lambda }_{n_s^u}]}^T$ and the semivariance function r, apart from using (7) to estimate the RSS value at position $(x_0,y_0)$, we can also estimate the uncertainty $S_0$, which is as follows:  

$$S_0=\sum _{n=1}^{n_s^u}{\lambda }_{n}r\left(d_{n0}\right),$$

(13)

where $d_{n0}$ is the distance between the sampled position $(x_n,y_n)$ and the unsampled position $(x_0,y_0)$. As shown in Figure 3, Figure 3a is the actual REM of the current area for UAV u at the stage s, Figure 3b represents the sampled points $n_s^u$ collected by other UAVs, Figure 3c shows the constructed REM based on the sampled points $n_s^u$, and Figure 3d depicts the uncertainty of each position.

In order to select appropriate sampling positions for UAV u at stage s, we normalize the RSS and uncertainty, respectively, and integrate them by adding with a weight $\eta$ to generate a probability map of sampling point distribution. The map indicates the probability distribution that UAV u will choose that positions for sampling, as shown in Figure 4. Since randomly located sampling points are suitable for parameter estimation [43], when UAV u adopts a granular sampling strategy at stage s, it will randomly select $(n_{set}-n_s^u)$ sampling positions based on the probability map. If $n_s^u$ is greater than $n_{set}$, UAV u does not perform sampling but instead proceeds directly to the next stage position $O_{s+1}$.

4.3. The Path Planning Method Based on DRL

After determining the starting position $O_s$, the next stage position $O_{s+1}$, and sampling positions for UAV u at stage s, a DRL method is proposed for planning the shortest path at stage s. We formulate the Markov decision process as follows:

• $State$: The $s_t$ represents a solution to the path planning at time step t, i.e., $s_t$= ($O_s$, $n_t^1$, ..., $n_t^i$, ..., $O_{s+1}$), where $n_t^i$ denotes the i-th sampling position for UAV u at stage s.
• $Action$: The action $a_t$ is represented by a point pair ($n_t^i$, $n_t^j$). We update the state $s_t$ using the 2-opt method [44].
• $Reward$: To optimize the initial solution as much as possible within time step T, we set the reward as follows:

  $$r_t=r(s_t,a_t,s_{s+1})=D\left(s_t^{*}\right)-min\{D\left(s_t^{*}\right),D\left(s_{t+1}\right)\},$$

  (14)

  where D(s) represents the sum of Euclidean distances between consecutive points in sequence s. $s_t^{*}$ is the best solution found before step t, which is updated only when $D\left(s_t^{*}\right)$ < $D\left(s_{t+1}\right)$. The cumulative reward for time step T to maximize is expressed as $G_T={\sum }_{t=0}^{T-1}{\gamma }^t{r}_t$, where $\gamma$ is the discount factor.

To train a policy ${\pi }_{\theta }$ to consistently generate reasonable actions $a_t$ at each time step, we employ the actor–critic algorithm. Inspired by [33], our actor network design is shown in Figure 5. Since positions $O_s$ and $O_{s+1}$ are fixed, the input to the actor network is the position $f\left(n_i\right)$ of each collecting point in the current state s. Since linear transformation alone cannot capture the position of each collecting position, we apply sinusoidal positional encoding to add the node position embedding information. Then, we use the Transformer [45] architecture to encode input features, obtaining the embedding result $h^1\left(n_i\right)$, with a dimensionality of 128. In order to effectively integrate global information into each node, we aggregate ${\left\{h^1\left(n_i\right)\right\}}_{i=1}^I$ by max-pooling to obtain $h_g$, embedding $h^1\left(n_i\right)$ into $h^c\left(n_i\right)$ through $h^c\left(n_i\right)$ = $W^1{h}^1\left(n_i\right)+W^g{h}_g$, where $W^1$, $W^g\in {\mathbb{R}}^{128\times 128}$. Given node embeddings $H^c=[h^c\left(n_1\right),\cdots ,h^c\left(n_I\right)]$, we obtain a compatibility matrix $Y\in {\mathbb{R}}^{I\times I}$ through a compatibility layer [45], which reflects the scores for selecting each pair of nodes. We transform Y into selection probabilities $P\in {\mathbb{R}}^{I\times I}$ for each node pair using a softmax layer. Since selecting the same two nodes is meaningless, we mask the diagonal elements. The elements $p_{ij}$ in matrix P represent the probability of selecting ($n^i$, $n^j$) for local operations.

The design of the critic network is similar to that of the actor network, with the following differences: (1) $h^c\left(n_i\right)$ is obtained through an average pooling layer, and (2) the score $v_{\varphi }$ is outputted through a fully connected layer. The complete algorithm is given in Algorithm 2.

Algorithm 2 Actor–Critic Algorithm-Based Path Planning

1: Input: Initial actor network parameters $\theta$, critic network parameters $\varphi$, number of epochs E, batches B, points M, step limit T, update step N. 2: Output: Trained policy ${\pi }_{\theta }$. 3: for each episode do 4:  randomly generate B batches of problem instances for path planning with M positions. 5:  for $b=1,\cdots ,B$ do 6:   select the problem instances for the b-th batch and randomly generate state $s_0$; 7:   $t\leftarrow 0$; 8:   while $t<T$ do 9:    Reset the gradients of $\theta$ and $\varphi$, $t_s=t$, and obtain state $s_t$. 10:    while $t-t_s<Nandt\ne T$ do 11:     sample $a_t$ based on ${\pi }_{\theta }\left(a_t\right|s_t)$; 12:     obtain reward $r_t$ and the next state $s_{t+1}$. 13:     $t\leftarrow t-1$ 14:    end while 15:    R = $v_{\varphi }\left(s_t\right)$; 16:    for $i\in \{t-1,\cdots ,t_s\}$ do 17:     $R\leftarrow r_i+\gamma R$; $\delta \leftarrow R-v_{\varphi }\left(s_i\right)$; 18:     $d\theta \leftarrow d\theta +{\sum }_{S_b}\delta \nabla log{\pi }_{\theta }\left(a_i\right|s_i)$; 19:     $d\varphi \leftarrow d\varphi +{\sum }_{S_b}\delta \nabla log{v}_{\varphi }\left(s_i\right)$ 20:    end for 21:    update $\theta$ and $\varphi$; 22:   end while 23:  end for 24: end for

In the actual deployment of the algorithm, since we directly use the trained policy to output the UAV’s flight path, the computational overhead mainly includes the forward propagation process. First, during the node embedding phase, we perform a linear transformation on the features of each node, which has a complexity of $\mathcal{O}(N·d·h_0)$, where N is the number of sampling positions for UAV u at stage s, d is the feature dimension, and $h_0$ is the embedding dimension. Additionally, the complexity of adding position encoding is $\mathcal{O}(N·h_0)$. Next, in the self-attention mechanism, we need to compute the relationships between nodes, resulting in a complexity of $\mathcal{O}(N^2·h_0)$, as we generate an $N\times N$ attention weight matrix and perform weighted summation. Subsequently, the complexity of action selection is $\mathcal{O}\left(1\right)$. By aggregating these analyses, we conclude that the overall computational complexity of the policy network is $\mathcal{O}(N^2·h_0)$. This complexity indicates that as the number of sampling positions for UAV u at stage s increases, the computational cost of the self-attention mechanism will significantly impact the efficiency of the algorithm. Therefore, the goal of our method is to construct a high-quality REM with a limited amount of RSS data. Subsequent experimental setups will also focus on discussing the quality of the REM constructed under low sampling rates.

5. Simulation Result

In this section, the performance of the proposed multi-stage approach to data sampling and path planning with multiple UAVs is evaluated. We compare the performance of our approach with other shortest path planning algorithms while also comparing the quality of the constructed REMs with other sampling strategies.

5.1. Simulation Setup and Parameter Settings

Simulation analysis and validation are conducted in an urban environment, where the RSS distributions are generated by the COST 231 extended WalfischIkegami empirical propagation model. The radiation source operates at a frequency of 2.4 GHz, with a transmission power of 43 dBm. Additionally, the target area is a square region with an area of 1000 × 1000 m², where the grid side length l is set to 5m, dividing the area into 200 × 200 small grids. Given the extensive area, three radiation sources are strategically positioned, each at an elevation of 30 m. The UAV is fixed at a height of H = 100 m to sample RSS data. When the sampling rate is too low, it is difficult to ensure the accuracy of REM construction, while a high sampling rate leads to excessively long sampling times. To demonstrate the superiority of our method in constructing high-quality REMs with minimal data, similar to other works on building REMs using UAV-sampled data, such as [24,40], the sampling rate $\epsilon$ for our experiments ranges from 0.25% to 2%, and $n_s^u$ is set to 12.5% of $\epsilon D^2$. The algorithm parameters are as shown in

Table 2. Algorithm parameters.

Parameter	Value	Parameter	Value
$R_{th}$	5 dB	$\alpha$	0.8
r	125 m	$\eta$	0.3

. In addition to experimenting with the REM generated through simulation, we also test the algorithm’s performance on two other datasets. The first dataset is the DeepREM dataset [46], which contains several urban scenario REMs. The other dataset is the Gudmundson dataset [8], which generates the RSS values using a path loss model and a shadowing model.

We used Pytorch 2.1.0 to implement our proposed algorithm, and all the codes are run on Microsoft Windows 11 with 13th Gen Intel(R) Core(TM) i9-13900KF CPU @3.00 GHz and NVIDIA GeForce RTX 4090 GPU. We trained 100 epochs for the policy $\pi$ of solving the shortest path between 50 points, with an initial learning rate of ${10}^{-4}$, a step limit T of 1000, and a batch size of 5120. After training the policy, we conducted testing with the policy’s step limit T set to 100.

5.2. Illustrative UAV Trajectories

The goal of path planning in this paper is to determine the shortest path for UAVs to traverse specified positions for accomplishing data sampling tasks. In Figure 6, we demonstrate the flight trajectories of the UAVs when U = 3 and U = 4, revealing a clear coordination and exploration among the UAVs. Each UAV departs from the starting position and is responsible for different parts of the target area during the data sampling mission. When the entire area is marked as explored, the UAVs will return. In Figure 6a, when U = 3, each UAV independently completes the data sampling of the region with dramatic RSS changes around the radiation source. However, in Figure 6b, when U = 4, the UAVs collaboratively perform the primary data sampling tasks in areas with dramatic RSS changes. These effects are a result of stage position selection based on Voronoi partitioning and a granular sampling strategy for choosing sampling positions.

5.3. Sampling Strategy Performance

Three construction methods including different sampling positions and different recovery methods were conducted via different sampling rates. The reconstruction method of this paper is denoted as Proposed-KGA, which indicates that after optimizing the sampling positions using the sampling strategy proposed in this paper, the KGA method is employed to reconstruct REM. The IDW algorithm [15] and the GPR algorithm [16] are also included as the representation of data-driven methods.

In order to compare the REM construction performance of the data sampled by our method, we compared it with the following sampling strategies:

Random sampling strategy [47]: Data are sampled at random positions, modeled by a homogeneous Poisson point process.

Grid sampling strategy [47]: Data are sampled at grid positions separated by a certain distance.

In the simulation environment, due to the sparse presence of buildings, the overall path loss is predominantly governed by free space loss. Regions where signal strength varies significantly are primarily located near the radiation sources. Thus, the distribution of sampling points in Figure 7b is mainly concentrated around the three radiation sources. In the Gudmundson dataset, the overall path loss incorporates a spatially correlated shadow as shadow fading on top of free space loss. Although the distribution of sampling points in Figure 7g is also around the radiation sources, it is not as densely distributed as in the simulation environment, but rather more scattered to capture the shadow fading effects. In the DeepREM dataset, which features a more complex urban environment with numerous buildings, significant variations in signal strength often result from building obstructions. Consequently, the distribution of sampling points in Figure 7l is scattered throughout the entire environment.

From Figure 7, it can be seen that the REM construction using the sampling strategy proposed in this paper exhibits better construction performance compared to other sampling strategies. This is because the proposed sampling strategy achieves a coordinated exploration of the entire area and increases the distribution of sampling positions in regions with dramatic RSS changes. The sampling strategy not only effectively restores the variations in RSS near the radiation source but also captures sudden changes in signal strength caused by obstacles, thereby enhancing the overall construction effectiveness.

To evaluate the impact of sampling strategies on the construction quality of the REM, this paper uses the following two metrics for assessment:

• Structural Similarity (SSIM): SSIM focuses on the structure and visual quality of the constructed REM, specifically whether it visually resembles the actual REM.
• Root Mean Square Error (RMSE): RMSE is defined as the difference between the actual REM values and the constructed REM values, given by

  $$RMSE=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{(Y_i-{\widehat{Y}}_i)}^2},$$

  (15)

  where $Y_i$ and ${\widehat{Y}}_i$ are the true and the estimated RSS values in dBm at the ith grid, respectively.

As shown in Figure 8, under the condition of a low sampling rate, the SSIM of the REM constructed by the method in this paper is higher than that of reconstruction methods based on other sampling strategies. This is attributed to the granular sampling strategy employed in regions with significant variations in RSS values, effectively capturing the large fluctuations in RSS caused by obstructions or proximity to the radiation source.

Intuitively, the larger the spatial sample density is, the more likely it is to find the data points that are close and hence highly correlated with the channel at the target location [47]. When the number of sampling points is fixed, the grid sampling strategy provides a higher spatial sampling density compared to other sampling strategies, resulting in the optimal RMSE when constructing the REM using interpolation methods. The sampling strategy proposed in this paper increases the spatial sampling density in regions where the RSS values change significantly due to obstruction by obstacles or proximity to radiation sources, making the RMSE of the REM constructed with this strategy close to that of the REM constructed with grid collecting strategy. The experimental results are reflected in Figure 9.

5.4. Path Planning Performance

The goal of path planning in this paper is to determine the shortest path for UAVs to traverse specified positions for accomplishing data sampling tasks. Regarding the DRL action $a_t$, we attempted 2-opt and node swap (i.e., directly swapping the positions of two nodes), and the rewards during the training process are shown in Figure 10. It can be observed that the final reward value obtained by 2-opt is higher.

We analyze the influence of each component in our actor network on the performance through ablation studies. The results are as shown in

Table 3. Ablation studies.

Self-Attention Layer	Compatibility Layer	Positional Encoding	Path Length	Reward
✓	✓	✓	6.137 ± 0.042	20.78
—	✓	✓	8.869 ± 0.086	18.51
✓	—	✓	9.293 ± 0.037	18.87
✓	✓	—	26.59  ± 0.005	9.16

Note: The ✓ indicates the inclusion of that specific component in the Actor Network.

. It can be seen that each component affects the performance of path planning.

To compare algorithm performance, we also implemented the following algorithms:

• Greedy algorithm: The greedy algorithm selects the nearest unsampled point from the current position of the UAV greedily, until all positions have been collected.
• Genetic algorithm (GA): It initializes several paths, with each path representing a chromosome. For each generation, it implements a crossover and mutation process. After several generations, the path with the shortest distance is considered the final solution.
• Simulated annealing (SA) algorithm: It sets an initial temperature for starting annealing and initializes a path as the initial solution. During each cooling process, it perturbs the current path to obtain a new path and updates the shortest path as the optimal solution. When the temperature decreases to the termination temperature level, it outputs the optimal solution.
• Recurrent neural network (RNN) [34]: The algorithm is based on a DRL framework, where a recurrent neural network is used to learn a construction heuristic policy, where each output is a position, gradually constructing a path.
• Attention model (AM) [35]: The algorithm is also used to learn a construction heuristic policy based on a DRL framework, and the network is an attention model.

For a fair comparison, multiple experiments are conducted to find the most suitable hyperparameters for the comparative algorithms. For the GA, the population size is set to 100, and we run for 50 generations. As for the SA, the initial temperature is set to ${10}^6$, with 500 iterations per cooling step and a cooling rate set to 0.98. Due to randomness, we independently conduct five rounds of experiments for each algorithm and average the output results. For the two DRL algorithms, RNN and AM, we start with the hyperparameters suggested in the original paper and continuously adjust them to train the highest-quality policy.

From Figure 11, it can be observed that as the sampling rate increases, our algorithm can find shorter flight paths compared to other methods. Figure 12 illustrates that the time required for running DRL algorithms is much shorter than for traditional evolutionary algorithms. This is because DRL algorithms use pre-trained policies during deployment, and the computation only involves forward propagation through the network. This indicates that our algorithm can find a shorter flight path within a lower time cost compared to other algorithms.

To demonstrate the scalability of our method, Figure 13 shows the trend of the average flight path length per UAV as the number of UAVs increases, with a fixed sampling rate of 1%. Due to the multi-UAV collaboration implemented by the method in this paper, as the number of UAVs increases, the number of positions each UAV needs to sample decreases, which in turn reduces the corresponding flight path length.

5.5. Impact of Parameter Settings

The parameters of our method mainly include the parameter $\alpha$ in (3), the threshold $R_{th}$, and the weight $\eta$ for generating the probability map.

For $\alpha$, different choices of $\alpha$ mainly affect the selection of stage positions. However, regardless of the value of $\alpha$, the entire area will always be explored, so its primary impact is on the length of the UAV’s flight path. The selection of $\alpha$ also needs to be designed differently in different environments. Figure 14 demonstrates the effect of $\alpha$ on the UAV’s flight path length in two environments: our simulation environment with fewer obstacles and the complex urban environment of the DeepREM dataset. In our simulation environment, as the value of $\alpha$ increases, the flight path length of the UAV decreases. This is because the area of signal strength variation is mainly near the radiation source, requiring the UAV to explore deeper into that region. However, in the complex urban environment, due to the impact of obstacles on signal strength, the UAV needs to strike a balance between deep exploration and collaborative exploration to avoid redundant exploration of the same areas.

To select an appropriate value for $R_{th}$, the SSIM values of the REMs constructed in two different environments are shown in Figure 15. As can be seen, in our simulation environment with fewer obstacles, a relatively larger $R_{th}$ value results in the highest SSIM. This is because the signal strength variations are primarily concentrated near the radiation source, and concentrating sampling points in these regions helps improve the accuracy of the REM construction. However, in the complex urban environment of the DeepREM dataset, a smaller $R_{th}$ value allows the UAV to focus on areas with significant signal variation caused by obstacles. However, a value that is too small leads to the UAV performing additional sampling in every region, which results in no available sampling points in later stages, thus affecting the overall reconstruction accuracy.

The purpose of the weight $\eta$ is to strike a balance between focusing on regions with stronger signal intensity and reducing uncertainty across the entire area. Figure 16 and Figure 17 illustrate the impact of different $\eta$ values on the accuracy of REM construction.

As shown in Figure 16, with an increase in the weight $\eta$, the SSIM value gradually increases. This is because larger $\eta$ values make the sampling strategy more focused on reducing uncertainty across the entire region, thereby improving the accuracy of REM construction. Figure 17 indicates that when the weight $\eta$ = 0.3, the RMSE reaches its minimum value. Sampling more points closer to the radiation source can effectively reduce the RMSE of REM construction [48]. Intuitively, regions with stronger signal intensity are usually closer to the radiation source. Therefore, appropriately increasing the number of sampling points in these areas helps reduce the RMSE of REM construction, thereby improving reconstruction accuracy.

6. Conclusions

In this paper, we propose a multi-stage approach to data sampling and path planning for multiple-UAV-enabled aerial REM construction. The collaboration among multiple UAVs is achieved by using Voronoi partitioning method. To improve the quality of constructing aerial REM, a granular sampling strategy is proposed, where the selection of sampling positions is based on the construction results and uncertainties. A self-attention-based policy network is designed and trained through the actor–critic framework to obtain a heuristic strategy for solving the shortest path problem for selected sampling positions. The simulation results illustrate that the algorithm achieves collaboration among multiple UAVs, reducing the length of UAV flight paths while improving the SSIM of the aerial REM and reducing the RMSE.