1030 Research

# Time-Varying Fuzzy Contour Trees

Anna-Pia Lohfink
e-mail<EMAIL_ADDRESS>Frederike Gartzky e-mail<EMAIL_ADDRESS>Florian Wetzels
e-mail<EMAIL_ADDRESS>Luisa Vollmer  Christoph Garth
e-mail<EMAIL_ADDRESS>Scientific Visualization Lab – Technische
Universität Kaiserslautern

###### Abstract

We present a holistic, topology-based visualization technique for spatial time
series data based on an adaptation of Fuzzy Contour Trees. Common analysis
approaches for time dependent scalar fields identify and track specific
features. To give a more general overview of the data, we extend Fuzzy Contour
Trees, from the visualization and simultaneous analysis of the topology of
multiple scalar fields, to time dependent scalar fields. The resulting time-
varying Fuzzy Contour Trees allow the comparison of multiple time steps that
are not required to be consecutive. We provide specific interaction and
navigation possibilities that allow the exploration of individual time steps
and time windows in addition to the behavior of the contour trees over all
time steps. To achieve this, we reduce an existing alignment to multiple sub-
alignments and adapt the Fuzzy Contour Tree-layout to continuously reflect
changes and similarities in the sub-alignments. We apply time-varying Fuzzy
Contour Trees to different real-world data sets and demonstrate their
usefulness.

Human-centered computingVisualizationVisualization techniquesTreemaps; Human-
centered computingVisualizationVisualization design and evaluation methods

Introduction

Analyzing time dependent scalar fields is usually done feature based by
identifying regions with specific patterns and tracking them between time
steps. While these approaches are helpful for specialized tasks, they are not
suited to provide a general overview of the complete data set due to their
focus on specific features and the chronology. We provide a holistic view of
the data set and its feature candidates by applying Fuzzy Contour Trees (FCTs,
[7]), giving a simultaneous overview of the topological structure of the
considered time steps. By defining sub-alignments that originate from the
_overall alignment_ containing all time steps, the development of the field
topology can be tracked, while also being able to consider the overall
topological structure of the data set and compare the structures of arbitrary
subsets, independent of the chronology.

After the introduction of FCTs as a tool for the simultaneous visualization of
the topological structure of ensembles, we present time-varying FCTs (T-FCTs),
specialized to the analysis of time series data. Since time series data
represents a series of individual results, the application of FCTs to this
problem is straight forward – instead of different ensemble members, we
consider different time points. However, application specific tasks require
further development of possible user interactions and visualization to ensure
usefulness and consistency. After discussing related work and providing
background on our previous work in Section 1, we make the following
contributions: We present the T-FCT interface in Section 2. Our new
contributions in the FCT back-end and front-end are discussed in Sections 3
and 4 respectively, with a short overview of the new interaction possibilities
in 5. We apply T-FCTs to different real-world data sets in Section 6 and
conclude with a short discussion on future work in Section 7.

##  1 Fuzzy Contour Trees and Related Work

Figure 1: The FCT framework developed by Lohfink et al.

We give an overview of previous work related with the novel contributions of
T-FCTs. Work that is related to FCTs in general, concerning contour trees,
their application in ensemble visualization and merging of graph-based
topological descriptors can be found in [7].

#### Time Dependent Contour Trees

A theoretical consideration of time-varying reeb graphs for continuous space-
time data was given by Edelsbrunner et al. [4]. Szymczak describes sub-domain
aware contour trees and uses them to track accumulated topological changes
between slices of the data set. While the evolution of iso-surfaces is
plotted, the contour trees are not visualized [13].

Also different interactive tools for the analysis of time varying contour
trees and iso-contours have been developed: An interactive exploration tool
for split/merge and contour trees for different time steps was developed by
Sohn et al. [12]. They define a topology change graph and use it to navigate
trees of individual time steps. Bajaj et al. provide multiple calculated
signature graphs on time-varying scalar fields. Interaction possibilities and
the additional visualization of single contour trees provide real-time exact
quantification in the visualization of iso-contours [1]. Kettner et al. take
this idea further to non-decomposable topological properties and higher
dimensions in the Safari interface [5]. Lukasczyk et al. define and visualize
spatio-temporal Reeb graphs in [8] to extract and visualize trajectories and
relationships of hotspots. Oesterling et al. show the evolution of extrema in
high-dimensional data by plotting a 1D landscape profile for each time step
and connecting peaks to indicate critical events. These events are determined
as structural changes in time-varying merge trees.

Other than T-FCTs, none of these approaches visualize contour trees for more
than a single time step or show the evolution of contour trees.

#### Feature Tracking in Time Dependent Data

Features have been identified and tracked in scalar fields in many different
ways. The components of interests can be superlevel sets [8, 12] or sub-
domains with special geometric and topological properties [2, 6, 11]. Multiple
features on different levels of interest are often tracked in tracking graphs
[8, 12, 9, 15].

In all of these cases, determined features are tracked between subsequent time
steps, tracing these features over time. In contrast, our approach provides an
overview over topological features instead of tracking individual features.
Considering the topological structure at all time steps and matching it, we
provide a more holistic view on the data. Similarities and differences of time
steps can be determined flexibly and independent of the adjacency of the time
steps.

#### Fuzzy Contour Trees

The foundation for T-FCTs was developed by Lohfink et al. with the Fuzzy
Contour Tree framework [7]. See Figure 1 for an overview. This framework
allows the simultaneous analysis of the topological structure of all members
of an ensemble at once. As underlying back-end functionality, the member
contour trees are matched using the tree alignment algorithm. This is done
incrementally: after aligning two contour trees, the resulting alignment and a
further contour tree are given as input for the alignment algorithm and so on.
Since the order of the contour trees has an impact on the result, the contour
trees are shuffled randomly before the execution.

The resulting alignment is a super tree of all contour trees. Thus, a layout
for the alignment can be transferred to the individual contour trees, yielding
common layouts for all of them. Superimposing the contour trees with this
layout results in a Fuzzy Contour Tree in its grouped layout. Bundling edges
at saddles and linking the edges’ opacity to the fraction of contour trees
that contain this specific edge results in the bundled layout. To further
improve the clarity of the visualization, Lohfink et al. implemented an
optimized branch spacing method that improves the positioning of branches on
their parent branch while keeping certain properties of the original FCT.
Multiple interaction possibilities are available that link the FCT to an
overview of the contained ensemble members. Both, FCTs and T-FCTs were
designed to support certain domain-agnostic, elementary analysis tasks. Due to
space limitations we reserve their detailed description and evaluation for
future work.

##  2 Interface Overview

Figure 2: The time-varying Fuzzy Contour Tree interface: Time selector and the
FCT of the selected time steps.

Similar to the original FCT described by Lohfink et al. in [7] that is applied
to ensembles, the T-FCT gives insight in the topological structure of time
dependent data by simultaneous visualization of multiple contour trees. In
addition, T-FCTs are tailored to the application to time dependent scalar
fields, considering different time steps of the same simulation or
measurement. This specialized application gives us the possibility to adapt
visualization of and interaction with the FCT to support the exploration of
time dependent fields at the best. Figure 2 shows the T-FCT interface. Both,
the time selector on top and the FCT below are enhanced versions of interface
components of the original FCT interface, providing specialized interaction
with time dependent data and ensuring a consistent visualization.

In the following we describe the background of this visualization and how to
generate it. Details on the T-FCT layout and the time selector can be found in
Sections 4 and 5 respectively. The T-FCT interface was implemented using
$D^{3}$, python and jupyter notebook.

##  3 The Time-Varying Alignment

Considering time dependent scalar fields, the main focus lies on a consistent
alignment over time, prohibiting a randomized order of the input contour
trees. Instead, the contour trees are aligned sequentially. Since the
structure of the first aligned time step imposes its structure on the whole
alignment, other choices for this time step might be beneficial depending on
the data set. Furthermore, the matching of nodes in consecutive time steps is
prioritized as follows: the scalar values of nodes in the alignment are set to
the value of the nodes in the lastly aligned contour tree, if they exist in
this contour tree. Else the value remains untouched. Thus, this contour tree
is effectively contained in the alignment with its scalar values. Hence,
aligning consecutive time steps behaves similar to matching adjacent trees,
enforcing the consistent matching of nodes. Since for the final layout, only
values from the individual time steps are used, changing the assigned values
in the alignment has no impact on node ordering.

The alignment process is very flexible in its application due to use of
different metrics. We implemented an overlap metric in addition to the
existing volume, persistence and combined metrics. It is defined as $1-J(A,B)$
where the Jaccard index $J(A,B)=\frac{|A\cap B|}{|A\cup B|}.$ The contour tree
alignment algorithm with optional matching over time is available in the TTK
development branch [14].

Analysis of time dependent data often takes place on the level of individual
time steps instead of the whole data set. Comparing the topological structure
of arbitrary sub-sets is rendered possible by the calculation of the
corresponding _sub-alignment_.

Sub-alignments are calculated based on an existing overall alignment. This
ensures consistent matching of nodes in different sub-alignments, avoids
problems with node identification between sub-alignments and guarantees a
consistent layout. To obtain a sub-alignment, contour trees of time steps that
are not selected are one by one subtracted from the alignment by decrementing
the frequency of all contained nodes. Then, nodes with frequency 0 are deleted
and connectivity of the result is restored. To do so, we determine parent -
child relations of the neighboring nodes by determining paths from these
neighbors to the fixed alignment root (see [7] for details) and connect the
nodes accordingly.

The resulting sub-alignments are alignments of the chosen sub-set of contour
trees. They are likely to be less optimal with respect to the chosen metric
than the result of the alignment algorithm heuristic would be. However, all we
need sub-alignments to be is a sub-tree of the overall alignment to have a
consistent and traceable layout. Employing the resulting FCT, analysis tasks
can be approached for the selected sub-set of time steps.

##  4 The Time-Varying Fuzzy Contour Tree

Figure 3: The trickle-down-layout of the T-FCT: The layout of the overall
alignment (Level I) carries over to sub-alignments (Level II) that are
compared in the time selector. Their layout propagates further to the
individual contour trees (Level III).

As a basis for the T-FCTs, the layout of the overall alignment is computed as
described in [7]. In the setting of T-FCTs, the cost function of the simulated
annealing doesn’t only take into account scalar values of leaves and the
persistence of branches, but also their existence in time. Only contemporary
branches (i.e. branches whose existence intervals overlap) with overlapping
bounding boxes are regarded as overlapping.

Similar to the layout of individual contour trees being based on the layout of
their alignment, the sub-alignments’ layout is based on the layout of the
overall alignment. An illustration for this trickle-down layout is given in
Figure 3.

To obtain the layout of a sub-alignment from the overall alignment, several
steps are taken: First, its branch decomposition is computed. Then, the order
of its branches is set according to the horizontal order in the overall
alignment. The main branch, as the only branch containing two leaf nodes, is
positioned according to its leaf node that can vary over time (the root node
always remains the same). This facilitates spotting changes in the main
branch.

With the horizontal order of the branches set, it is optionally adapted to
avoid gaps that occur if branches from the overall alignment are not present
in the sub-alignment. While this results in more uniform individual layouts,
tracking the development of individual branches is easier without re-
positioning, especially in animations.

After the layout for the sub-alignment is determined, leaves, saddles and
branches are transferred from the individual contour trees as described in
[7], following paths in the sub-alignment. The coloring of leaf nodes is kept
consistent over different sub-alignments.

In contrast to common feature tracking methods, T-FCTs allow the direct
comparison of non-consecutive time steps as well as a topological overview
over all time steps. Matching order and metric are adaptable, providing in-
depth control of the result.

##  5 Interaction

Due to space limitations, we present added and enhanced interaction
possibilities only shortly. For more details see the FCT paper [7] and the
supplemental material.

Figure 4: The time selector: The complete data set with highlighted selected
time steps. Buttons allow the manipulation of the selection.

Interaction possibilities with the FCT aim on linking information concerning
individual contour trees with information on the overall behavior of the
individual fields. To adapt FCTs to the time-varying setting, the grid
providing information on individual contour trees and enabling tree and member
highlighting was turned into the time selector (Figure 4, compare Figure 1).
Instead of showing only time steps that are contained in the current sub-
alignment, every time step is represented by a colored slice of the selector.
Selected time steps are highlighted as boxes with the number of the time step.

To indicate areas of interest, the time selector is colored based on different
measures on the sub-alignment containing a time window of given width with the
considered time step at the center. Either these values are shown directly, or
the averaged sum of their node-wise differences with the neighboring time step
to illustrate changes between them. We implemented two centrality measures:
degree and betweenness centrality [10], focusing on edge insertions and
deletions, and structural changes respectively. Examples for the different
coloring options are given in Figure 5.

Figure 5: Time selector options Different selections and coloring.

Time steps can be selected via the time selector in three different modes:
_Window selection_ , _multi-selection_ and _periodic selection_ as special
case of multi-selection. See Figure 5 for examples.

_Window selection_ selects a connected sub-interval centered at the selected
time step. _Multi-selection_ allows the manual selection of multiple, not
necessarily adjacent time steps to provide deeper insight in individual
members and their differences. As a special case of multi-selection, the
_periodic selection_ automatically selects time steps with a given period and
shows a periodic marker.

A given selection can be shifted left and right. A play and a pause button
start and stop automated shifting to the right (Figure 4).

Figure 6: Tree highlighting: Right: Hovering selected time steps highlights
the individual contour tree in the FCT. Left: For time steps that are not part
of the current sub-alignment, only contained branches are highlighted without
indication of specific nodes and saddles.

_Tree highlighting_ is also generalized for time steps that are currently not
selected: Branches that are present are highlighted without marking specific
saddles and leaves (see Figure 6).

_Member highlighting_ in the time selector does not only take place for
currently selected time steps, but for all available time steps. Like this,
navigating the data set and finding patterns in branch occurrences is
facilitated. See Figure 7 for an example.

##  6 Applications

Figure 7: Periodic behavior in the sea ice data set. The marked branch is
contained in the contour tree only during winter months. (period 12, volume
metric)

###  6.1 Sea Ice

The sea ice data set describes arctic sea ice concentrations and is provided
by the National Snow and Ice Data Center [3]. The concentrations are given as
tenths of grid square area covered by ice. We set the value -1 for grid cells
that are not over sea. We applied T-FCTs to sea ice concentrations from 1980
to 1995 and applied a mild smoothing filter to the data to avoid non-binary
contour trees.

Containing only discrete values, large plateaus in the data result in contour
trees with branches of persistence 0, forbidding the use of metrics that rely
on persistence in the alignment process. Hence, we chose volume and overlap
metric for this example. The discrete nature of the saddles and extrema also
complicates finding a suitable layout for the T-FCT. Optimized branch spacing
as described in [7] shifts the branches to still allow a clear layout.

Figures 7 and 9 show the T-FCT for the sea ice data set with branches of
typical behavior highlighted. The selected branch in Figure 7 can be found in
the winter months of every selected year, reflecting the periodical increase
of sea ice during winter months. The highlighted minimum (ocean) is only
matched if it is surrounded by sufficiently large maxima (sea ice). In Figure
9 on the other hand, the selected branch vanishes after 1986, to re-occur in
early 1994. This branch represents an area of the arctic sea that was
frequently surrounded by sea ice during winters before around 1985 but not any
more. Both of these behaviors, periodical occurrence over the whole data set
and periodical occurrence only in about the first half of the considered time
frame with potential re-occurrence between 1992 and 1994 can be seen in
multiple branches.

Periodic selection allows the selection of individual months over all years.
An example for the comparison of individual months from 1980 to 84 is given in
Figure 8. Here, we used the overlap metric for the alignment and illustrate
the usefulness of optimized branch spacing in discrete data sets. The much
higher complexity of the main structures in the FCT during winter months
indicates the larger extent of sea ice in different ice floes. Although these
floes are interrupted by islands and potentially even free to move,
complicating a matching between time steps, similar structures are visible
between the months, indicating areas of the arctic ocean that are covered by
ice around the whole year.

Figure 8: Sea ice extent during individual months. Similar structures are
highlighted. For March, the comparison between the FCT with and without
optimized branch spacing is given (overlap metric). Figure 9: Vanishing
branches over time show the yearly decreasing extent of the sea ice. In 1994
the sea ice extent peaked and the branch re-occurs (period 12, volume metric).

###  6.2 Convection Simulation

Flow around a hot pole over time was simulated. The resulting 2D flow field is
of size 128 × 256. Material at rest is heated around the pole, begins to rise,
and forms a plume. An exemplary member is shown on the right of Figure 1. This
behavior is clearly reflected in the T-FCTs in Figure 10. Advancing the
selected window step by step through the data provides an overview of the
topological structure at individual time steps but also their connection and
development: Until time step 74, the pole is heating up. Then, the plume forms
and more and more temperature minima are enclosed by the formed plume. The
time selector color shows the degree values, reflecting the dynamic behavior
of the plume towards the end. We chose representative time steps and scaled
them to highlight the raising temperature; find the complete animation in the
supplemental material.

Figure 10: Sliding Window: Advancing step by step through the convection
simulation result provides a clear understanding of the ongoing processes. We
chose an expressive subset of time steps, the step number is given above the
FCTs.

##  7 Conclusion and Future Work

In this paper, we applied FCTs in the time-varying setting. As our new
contribution, we adapted back-end and front-end to its specific challenges,
resulting in the T-FCT interface. We illustrated usefulness and limitations in
real-world examples. As future work, different options to deal with cluttered
(Fuzzy) contour trees could be included, for example hiding nondescript
branches and detail on demand. Also, instead of a purely additive comparison
of contour trees, a subtractive comparison just showing branches not contained
everywhere could improve the overview. In summary, extending and applying FCTs
to time varying data results in useful, coherent visualizations providing a
more general overview of the data than common analysis approaches.

###### Acknowledgements.

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Re- search
Foundation) – 252408385 – IRTG 2057

## References

  * [1] C. L. Bajaj, V. Pascucci, and D. R. Schikore.  The contour spectrum.  In Proceedings. Visualization ’97 (Cat. No. 97CB36155), pp. 167–173, 1997. doi: 10 . 1109/VISUAL . 1997 . 663875
  * [2] P.-T. Bremer, E. M. Bringa, M. A. Duchaineau, A. G. Gyulassy, D. Laney, A. Mascarenhas, and V. Pascucci.  Topological feature extraction and tracking.  Journal of Physics: Conference Series, 78:012007, jul 2007. doi: 10 . 1088/1742-6596/78/1/012007
  * [3] W. L. Chapman and J. E. Walsh.  Arctic and southern ocean sea ice concentrations, version 1, arctic.  Boulder, Colorado USA. NSIDC: National Snow and Ice Data Center, 1991, updated 1996, Accessed: 2021-02-20. doi: 10 . 7265/N5057CVT
  * [4] H. Edelsbrunner, J. Harer, A. Mascarenhas, V. Pascucci, and J. Snoeyink.  Time-varying reeb graphs for continuous space–time data.  Computational Geometry, 41(3):149–166, 2008. doi: 10 . 1016/j . comgeo . 2007 . 11 . 001
  * [5] L. Kettner, J. Rossignac, and J. Snoeyink.  The safari interface for visualizing time-dependent volume data using iso-surfaces and contour spectra.  Computational Geometry, 25(1):97–116, 2003.  European Workshop on Computational Geometry - CG01. doi: 10 . 1016/S0925-7721(02)00132-3
  * [6] T. Lee, X. Tong, H. Shen, P. C. Wong, S. Hagos, and L. R. Leung.  Feature tracking and visualization of the madden-julian oscillation in climate simulation.  IEEE Computer Graphics and Applications, 33(4):29–37, 2013. doi: 10 . 1109/MCG . 2013 . 36
  * [7] A.-P. Lohfink, F. Wetzels, J. Lukasczyk, G. H. Weber, and C. Garth.  Fuzzy contour trees: Alignment and joint layout of multiple contour trees.  Computer Graphics Forum, 39(3):343–355, 2020. doi: 10 . 1111/cgf . 13985
  * [8] J. Lukasczyk, R. Maciejewski, C. Garth, and H. Hagen.  Understanding hotspots: A topological visual analytics approach.  In Proceedings of the 23rd SIGSPATIAL International Conference on Advances in Geographic Information Systems. Association for Computing Machinery, 2015. doi: 10 . 1145/2820783 . 2820817
  * [9] J. Lukasczyk, G. Weber, R. Maciejewski, C. Garth, and H. Leitte.  Nested tracking graphs.  Computer Graphics Forum, 36(3):12–22, 2017. doi: 10 . 1111/cgf . 13164
  * [10] M. Roy, S. Schmid, and G. Trédan.  Modeling and measuring graph similarity: The case for centrality distance.  In Proceedings of the 10th ACM International Workshop on Foundations of Mobile Computing, FOMC ’14, p. 47–52. Association for Computing Machinery, 2014. doi: 10 . 1145/2634274 . 2634277
  * [11] A. Schnorr, D. N. Helmrich, D. Denker, T. W. Kuhlen, and B. Hentschel.  Feature tracking by two-step optimization.  IEEE Transactions on Visualization and Computer Graphics, 26(6):2219–2233, 2020. doi: 10 . 1109/TVCG . 2018 . 2883630
  * [12] B. Sohn and C. L. Bajaj.  Time-varying contour topology.  IEEE Transactions on Visualization and Computer Graphics, 12(1):14–25, 2006. doi: 10 . 1109/TVCG . 2006 . 16
  * [13] A. Szymczak.  Subdomain aware contour trees and contour evolution in time-dependent scalar fields.  In International Conference on Shape Modeling and Applications 2005 (SMI’ 05), pp. 136–144, 2005. doi: 10 . 1109/SMI . 2005 . 45
  * [14] J. Tierny, G. Favelier, J. A. Levine, C. Gueunet, and M. Michaux.  The topology toolkit.  IEEE Transactions on Visualization and Computer Graphics, 24(1):832–842, 2018. doi: 10 . 1109/TVCG . 2017 . 2743938
  * [15] W. Widanagamaachchi, C. Christensen, V. Pascucci, and P. Bremer.  Interactive exploration of large-scale time-varying data using dynamic tracking graphs.  In IEEE Symposium on Large Data Analysis and Visualization (LDAV), pp. 9–17, 2012. doi: 10 . 1109/LDAV . 2012 . 6378962


# Trajectory and Power Design for Aerial CRNs with Colluding Eavesdroppers
††thanks: Manuscript received.

Hongjiang Lei, Jiacheng Jiang, Haosi Yang, Ki-Hong Park,
Imran Shafique Ansari, Gaofeng Pan, and Mohamed-Slim Alouini This work was
supported by the National Natural Science Foundation of China under Grant
61971080.H. Lei is with the School of Communications and Information
Engineering, Chongqing University of Posts and Telecommunications, Chongqing
400065, China and also is Chongqing Key Lab of Mobile Communications
Technology, Chongqing 400065, China (e-mail: leihj@cqupt.edu.cn).J. Jiang, and
H. Yang are with the School of Communications and Information Engineering,
Chongqing University of Posts and Telecommunications, Chongqing 400065, China
(e-mail<EMAIL_ADDRESS>cquptyhs@163.com).I. S. Ansari Bharath Institute of
Science and Technology, Chennai 600073, India (e-mail: ansarimran@ieee.org).G.
Pan is with the School of Cyberspace Science and Technology, Beijing Institute
of Technology, Beijing 100081, China (e-mail: gaofeng.pan.cn@ieee.org).K.-H.
Park and M.-S. Alouini are with CEMSE Division, King Abdullah University of
Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia (e-mail:
<EMAIL_ADDRESS>slim.alouini@kaust.edu.sa).

###### Abstract

Unmanned aerial vehicles (UAVs) can provide wireless access services to
terrestrial users without geographical limitations and will become an
essential part of the future communication system. However, the openness of
wireless channels and the mobility of UAVs make the security of UAV-based
communication systems particularly challenging. This work investigates the
security of aerial cognitive radio networks (CRNs) with multiple uncertainties
colluding eavesdroppers. A cognitive aerial base station transmits messages to
cognitive terrestrial users using the spectrum resource of the primary users.
All secondary terrestrial users and illegitimate receivers jointly decode the
received message. The average secrecy rate of the aerial CRNs is maximized by
jointly optimizing the UAV’s trajectory and transmission power. An iterative
algorithm based on block coordinate descent and successive convex
approximation is proposed to solve the non-convex mixed-variable optimization
problem. Numerical results verify the effectiveness of our proposed algorithm
and show that our scheme improves the secrecy performance of airborne CRNs.

###### Index Terms:

Unmanned aerial vehicle, trajectory and power design, average secrecy rate,
multiple colluding eavesdroppers.

##  I Introduction

###  I-A Background and Related Works

Due to their high mobility, low cost, and on-demand deployment flexibility,
unmanned aerial vehicles (UAVs) are considered critical to the future of the
Internet of Things (IoT), enabling the interconnection of all things [1, 2].
More and more communication scenarios are now becoming flexible and diverse
due to the inclusion of UAVs, such as UAV-assisted real-time surveillance,
traffic control, and emergency disaster relief. In these application
scenarios, UAV-assisted communication is considered to be an essential
component of future mobile communication systems. By designing the position of
UAVs or optimizing their flight trajectory, the air-to-ground (A2G) links
between UAVs and terrestrial nodes are line-of-sight (LoS) with a high
probability. Therefore, the trajectory of UAVs has a decisive impact on the
performance of the UAV-based communication system and becomes a fundamental
problem to be solved in the design of airborne communication systems [3].

The concept of the “Internet of Everything” has led to an increasingly
prominent role for UAVs in wireless communication systems [4]. UAV-assisted
communication systems utilize drones to provide ubiquitous wireless coverage
for a given area. In varying application scenarios, drones can be used in
different roles in communication systems, such as base stations (BSs) [5] \-
[8], relays [9] \- [11], and users [12] \- [14]. In [5], a communication
system consisting of multiple aerial BS was studied. The UAV trajectory was
optimized to maximize the minimum throughput under the consideration of user
scheduling coefficient and UAV transmit power. An efficient iterative
algorithm based on successive convex approximation (SCA) and block coordinate
descent (BCD) was proposed for mixed integer non-convex optimization. Their
results showed that the proposed algorithm could improve the system’s utility.
Considering the limited energy consumption of UAVs’ flight and ensuring the
system’s service quality, the author in [6] investigated the max-min fairness
problem to maximize the throughput of the system by jointly optimizing the UAV
transmit power and flight trajectory. In [7], a space-air-ground three-tier
heterogeneous network consisted of a satellite, multiple UAVs, and multiple
terrestrial BS to provide seamless coverage service for the ground users. A
two-stage joint optimization was proposed to deal with the intra- and inter-
interference and cross-tier resource allocation among different networks. A
multi-UAVs-aided mobile edge computing system was designed to balance traffic
and a new multi-UAV deployment mechanism was proposed in [8]. A near-optimal
solution algorithm was proposed to solve the load-balancing multi-UAV
deployment problem. When UAVs are used as relay nodes, long-distance wireless
links between terrestrial BSs and terrestrial nodes (TNs) can be quickly
established, which is helpful in disaster or other emergency scenarios [1]. In
[9], authors considered a UAV-assisted cooperative system in which multiple
UAVs were utilized as decode-and-forward relays to forward signals to the TNs.
The minimum transmission rate of the system was maximized by jointly
optimizing the flight trajectory of the relay UAVs, the transmit power on UAVs
and terrestrial BS. An iterative algorithm was proposed to solve the joint
optimization problem based on the block coordinate ascent technique,
introducing the slack variable and SCA techniques. In [10], an aerial relay
balanced the difference in transmission rates between free-space optical and
radio frequency links. The throughput for delay-limited and delay-tolerant
scenarios was maximized by carefully designing the flight trajectory of the
aerial relay, which in turn maximized the throughput. Combining intelligent
reflecting surface (IRS) and UAV can establish new channels or improve the A2G
quality through full-angle reflection in the air. A UAV-IRS-aided covert
scheme was proposed in [11], where the UAV carried an IRS and was utilized as
a relay. Considering Willie’s imperfect position state information, the covert
transmission rate was maximized by alternately optimizing Alice’s transmit
power, IRS phase shift, and the horizontal position of the UAV-IRS. In UAV-
aided wireless sensor networks and IoT communication systems, UAVs are
utilized as aerial access points to disseminate/collect information to/from
TNs [1]. A mixed-integer non-convex optimization problem was proposed to
prolong the wireless sensor network lifetime and the UAV’s trajectory was
jointly optimized to minimize the maximum energy consumption of all sensor
nodes in [12]. Refs. [13, 14] combined UAV trajectory design and wireless
charging technology in the scenarios of the Internet of Things, considered the
two-way link between information transmission and energy acquisition, and
realized the problem of maximum and minimum data acquisition by optimizing the
user scheduling coefficient and UAV flight trajectory.

The security of UAV communication applications is vulnerable to threats due to
the broadcast nature and the prevailing LoS channel conditions [15], [16],
[17]. Since then, more and more researches have focused on the security of
UAV-based communication systems. The security of A2G and ground-to-air (G2A)
links was investigated in [18], and the average secrecy rate (ASR) were
maximized by jointly optimizing the transmission power and the trajectory of
the UAV. By using successive convex approximation (SCA) and block coordinate
descent (BCD), iterative algorithms were proposed to solve the non-convex
problems. In [19], the security of underlay cognitive radio networks (CRNs)
was investigated. The ASR was maximized by designing an antenna jammer’s
trajectory and transmit power. However, the eavesdroppers were assumed to
operate in half-duplex mode and their position was assumed to be perfectly
known.

The full-duplex (FD) eavesdroppers pose more severe security threats to UAV-
based communication systems because they can not only intercept confidential
information from A2G links but also send jamming signals to degrade the
reception quality of the legitimate link. In [20], authors analyzed the
security risks caused by terrestrial FD eavesdroppers and propose an iterative
algorithm to maximize the ASRs of the downlink and uplink transmissions by
optimizing the UAV flight trajectory and transmission power. The proposed non-
convex problem was decomposed into sub-problems in BCD scheme with suboptimal
solutions obtained by the SCA. The position of the eavesdroppers at the base
stations was also assumed to be perfectly known.

Estimating the exact location and channel state information in these scenarios
is very challenging, especially for malicious eavesdroppers with some
deception equipment. In these scenarios, the bounded eavesdropper’s location
error model was utilized and the worst ASR (WASR) was considered. Ref. [21]
investigated the security of an A2G communication system with multiple
eavesdroppers, whose locations are assumed to be imperfect. Another UAV was
utilized as a friendly jammer to enhance the security of the UAV-aided
communication, Ref. [22] designed both the aerial BS and the aerial jammer’s
transmit power and trajectories to maximize the WASR. The WASR and the secrecy
energy efficiency (SEE) of the considered system were maximized by jointly
optimizing the trajectory and transmit power of the aerial base station,
respectively. Considering the imperfect location of the eavesdroppers and the
practical propulsion energy consumption, two cooperative dual UAV-assisted
data collection schemes have been proposed in [23]. The WASR and energy
efficiency were maximized by optimizing the flight trajectory of the dual
UAVs, the transmission power, and the scheduling coefficients. Furthermore, in
some practical scenarios, multiple eavesdroppers attempt to intercept
confidential information. In these scenarios, the base station needs to move
away from all the eavesdroppers as far as possible and reach the legitimate
targets as closely as possible to increase the secrecy rate. Authors in [24]
studied the secrecy performance of the energy-constrained dual-UAV systems
with multiple inaccurate location eavesdroppers. The WASR was maximized by
optimizing the 3D trajectories, the transmit power of the airborne transmitter
and jammer, and the time slot allocation scheme of the airborne transmitter.
The security of an underlay CRN with multiple imprecise location eavesdroppers
was investigated in [25]. The ASR of the cognitive user was maximized by
optimizing the UAV’s trajectory and transmission power. Authors in [26]
investigated the secure CRNs and a friendly jammer was utilized to transmit
artificial noise. Both scenarios were considered in which all the terrestrial
nodes’ locations were known and unknown. A dual collaborative underlay aerial
IoT system with multiple uncertainty potential eavesdropper was considered and
the WASR was maximized by jointly considering UAV trajectories, user
scheduling, and transmit power in [28].

To maximize the eavesdropped information in some scenarios, the eavesdroppers
may engage in colluding behavior, i.e., gathering and sharing information,
especially in those scenarios where the quality of the eavesdropping links is
poor. The most hazardous scenario is when all eavesdroppers are involved in
malicious collusion [27]. The security of cellular UAV communication systems
with multiple legitimate and illegitimate ground receivers was investigated in
[29]. The secrecy rate of the UAV communication systems was maximized for both
non-colluding and colluding scenarios by jointly optimizing the UAVs’ 3D
location and transmit power. However, the location of all eavesdroppers was
assumed to be perfectly known. Authors in [30] studied the energy efficiency
of a UAV-aided system with multiple colluding eavesdroppers with imperfect
location information and formulated a mixed-integer non-convex optimization
problem to minimize the energy consumption by jointly optimizing the UAV’s
trajectory and user scheduling. Authors in [31] investigated the secrecy
performance of a UAV-aided non-orthogonal multiple access (NOMA) systems with
multiple colluding eavesdroppers with imperfect location information. The
minimum WASR was maximized by jointly optimizing the UAV’s trajectory, the
transmit power, and the power allocation coefficient. Table I outlines recent
works in literature related to trajectory design of UAV communication systems.

TABLE I: Related to Trajectory Design of UAV Communication Systems. Reference | Security | Role of UAVs |  Imperfect Location of Eavesdroppers |  Colluding Eavesdroppers | FD Eavesdroppers | CRNs
---|---|---|---|---|---|---
[5],[6] |  | BS |  |  |  |
[9]-[11] |  | Relay |  |  |  |
[12]-[14] |  | User |  |  |  |
[18] | ✓ | BS, User |  |  |  |
[19] | ✓ | Jammer |  |  |  | ✓
[20] | ✓ | BS, User |  |  | ✓ |
[21] | ✓ | BS | ✓ |  |  |
[22] | ✓ | BS, Jammer | ✓ |  |  |
[23] | ✓ | User, Jammer | ✓ |  |  |
[24] | ✓ | BS, Jammer | ✓ |  |  |
[25] | ✓ | BS | ✓ |  |  | ✓
[26] | ✓ | Jammer | ✓ |  |  | ✓
[27] | ✓ | BS |  | ✓ |  |
[28] | ✓ | BS, Jammer | ✓ |  |  | ✓
[29] | ✓ | BS |  | ✓ | ✓ |
[30] | ✓ | BS | ✓ | ✓ |  |
[31] | ✓ | BS | ✓ | ✓ |  |
Our Work | ✓ | BS | ✓ | ✓ | ✓ | ✓

###  I-B Motivation and Contributions

The previous discussed works demonstrated that the security of UAV-aided CRNs
can be significantly improved by designing the UAV’s trajectory and other
system parameters (such as transmit power). However, these outstanding works
did not answer the following security issue:  How to design the trajectory and
transmit power of the aerial base station where multiple inaccurate location
FD eavesdroppers intercept information in colluding mode? What is the optimal
secrecy performance of the aerial CRNs with multiple colluded FD eavesdroppers
with imperfect location?  Hence, this work answers these questions by
optimizing the WASR of aerial CRNs. It is assumed that multiple location-
uncertainty illegitimate receivers work in FD and colluding modes and the
legitimate receivers jointly decode the received legitimate messages via
maximal ratio combining (MRC) scheme. The main contributions of this work are
summarized as follows:

  1. 1.

We consider an underlay CRN consisting of an aerial base station and multiple
primary and cognitive users. Multiple malicious FD eavesdroppers at uncertain
locations wiretap confidential information and transmit interference to
degrade the reception quality of the legitimate links. All terrestrial
secondary users and illegitimate receivers are assumed to utilize the
cooperative decoding of received message, respectively. The WASR is maximized
by jointly optimizing the UAV’s transmission power and trajectory while
satisfying SEE limit. Due to its non-convexity, the optimization problem is
split into two sub-problems and transformed into approximated convex forms via
SCA. The BCD technique is utilized to solve these subproblems.

  2. 2.

The proposed algorithm’s convergence and complexity are analyzed and the
numerical results of the proposed scheme are compared with three benchmarks
where the aerial base station works with the fixed trajectory and transmits
power simultaneously or only optimizing the trajectory and transmit power
separately. The efficiency and the convergence of the proposed algorithm are
verified.

  3. 3.

Relative to [19], [25], [26, 28] where the secure CRNs with a single
eavesdropper was considered, this work investigates the security of aerial
CRNs with multiple colluding eavesdroppers residing at uncertain location.
Technically speaking, considering the collusion scenarios is much more
challenging than the independent eavesdropping scenarios.

  4. 4.

Although the colluding scenarios were considered in [27] and [29], the results
can not be applied to CRNs directly since the location of the eavesdroppers is
assumed to be perfectly known and the underlay condition was not considered.
This work assumes that both legitimate and illegitimate receivers are
considered jointly to decode the received legitimate messages. Both the
uncertainty of eavesdroppers’ locations and the underlay condition are
considered simultaneously.

###  I-C Organization

The rest of this paper is organized as follows. The system model and problem
formulation are provided in Section II and III, respectively. A two-step
alternating algorithm is proposed in Section IV to solve the problem.
Simulation results are demonstrated in Section V. Finally, Section VI
concludes this paper. TABLE II lists the notations and symbols utilized in
this work.

TABLE II: List of Notations. Notation | Description
---|---
${\mathbf{q}}_{S}^{0}$ | Initial location of $S$
${\mathbf{q}}_{S}^{F}$ | Final location of $S$
${{\mathbf{w}}_{{D_{k}}}}$ | Horizontal location of $D_{k}$
${{\mathbf{w}}_{{E_{m}}}},{{\mathbf{\hat{w}}}_{{E_{m}}}}$ | Exact and estimated location of $E_{m}$
${{\mathbf{w}}_{{U_{r}}}}$ | Horizontal location of $U_{r}$
${r_{{E_{m}}}}$ | Maximized estimation errors of the distance
$H$ | Altitude of $S$
${\Gamma_{r}}$ | Tolerance threshold of $U_{r}$
${\sigma^{2}}$ | The noise power
$P_{S}^{\max}$ | The peak transmit power of $S$
${P_{{E}}}$ | Transmit power of $E_{m}$
$\alpha$ | Path loss exponent
$V_{S}^{\max}$ | Maximum flight speed of the UAV
${\beta_{0}}$ | Channel power gain at the reference distance
${\delta_{t}}$ | Time slot length
$\varepsilon$ | Algorithm convergence precision
$\Psi$ | SEE threshold of $S$

##  II System Model

Figure 1: UAV-enabled CRN consists of an aerial base station (${S}$), multiple
primary users (${U_{r}}$), and multiple cognitive users (${D_{k}}$). Multiple
FD eavesdroppers (${E_{m}}$) at uncertain locations not only try to wiretap
the confidential information but also send interference to degenerate the
reception quality of the legitimate link.

As shown in Fig. 1, we consider an aerial underlay communication system
consisting of a rotary-wing UAV that works as a base station (${S}$) and
multiple terrestrial users (${D_{k}},k=1,\cdots,K$) in the presence of
multiple primary users (${U_{r}},r=1,\cdots,R$) and multiple malicious FD
eavesdroppers (${E_{m}},m=1,\cdots,M$). All the nodes are equipped with a
single antenna.

To avoid collisions, $S$ flies at a constant altitude $H$ that is higher than
the highest obstacle in the service area [21], [23], [25], [32]. The flight
period, $T$, is divided into ${N}$ time slots as ${\delta_{t}}=\frac{T}{N}$
[3]. Without loss of generality, a 3D cartesian coordinate is used. Then the
coordinates of $S$, $D_{k}$, ${U_{r}}$, and $E_{m}$ are written as
${{\mathbf{q}}_{S}}\left(n\right)={\left[{{x_{S}}\left(n\right),{y_{S}}\left(n\right)}\right]^{T}}$,
${{\mathbf{w}}_{{D_{k}}}}={\left[{{x_{{D_{k}}}},{y_{{D_{k}}}}}\right]^{T}}$,
${{\mathbf{w}}_{{U_{r}}}}={\left[{{x_{{U_{r}}}},{y_{{U_{r}}}}}\right]^{T}}$,
and
${{\mathbf{w}}_{{E_{m}}}}{\left[{{x_{{E_{m}}}},{y_{{E_{m}}}}}\right]^{T}}$,
respectively, where $n=1,\cdots,N$.

Similar to [19] and [23], the A2G links are assumed to be LoS links and the
channel coefficient between $S$ and the terrestrial receiver $i$ is expressed
as

${h_{Si}}\left(n\right)=\frac{{{\beta_{0}}}}{{d_{Si}^{2}\left(n\right)+{H^{2}}}},$
(1)

where ${i\in\left\\{{{D_{k}},{U_{r}},{E_{m}}}\right\\}}$, ${\beta_{0}}$
denotes the channel power gain at the reference distance and
${d_{Si}}\left(n\right)=\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{i}}}\right\|$
signifies the Euclidean distance between $S$ and $i$, respectively.

Due to the uncertainty of the eavesdroppers’ locations, ${{h_{SE}}}$ obtained
by (1) is imperfect. We define
${{\mathbf{\hat{w}}}_{{E_{m}}}}\in{\mathbb{R}^{2\times 1}}$ and
${r_{{E_{m}}}}$ to express the estimated location and maximum error of the
difference between the estimation and the practical distance, respectively.
According to [23], it is assumed that
$\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{{\mathbf{\hat{w}}}}_{E}}}\right\|\geq{r_{E}}$.
Following the triangle inequality, ${d_{S{E_{m}}}}$ is expressed as

$\displaystyle{d_{S{E_{m}}}}\left(n\right)$
$\displaystyle=\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{{E_{m}}}}}\right\|$
(2)
$\displaystyle\geq\left|{\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{{\mathbf{\hat{w}}}}_{{E_{m}}}}}\right\|-\left\|{{{{\mathbf{\hat{w}}}}_{{E_{m}}}}-{{\mathbf{w}}_{{E_{m}}}}}\right\|}\right|$
$\displaystyle\geq\left|{\underbrace{\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{{\mathbf{\hat{w}}}}_{{E_{m}}}}}\right\|}_{\buildrel\Delta\over{=}{{\hat{d}}_{S{E_{m}}}}\left(n\right)}-{r_{{E_{m}}}}}\right|,$

where ${{\hat{d}}_{S{E_{m}}}}\left(n\right)$ denotes the estimation of
${d_{S{E_{m}}}}\left(n\right)$.

The ground-to-ground (G2G) link between $D_{k}$ and $E_{m}$ is assumed to
undergo quasi-static Rayleigh model, then channel coefficient is expressed as

${h_{{D_{k}}{E_{m}}}}=\frac{{{{\beta_{0}}}}}{{{d_{{D_{k}}{E_{m}}}^{\alpha}}}}{\xi_{{D_{k}}{E_{m}}}},$
(3)

where
${d_{{D_{k}}{E_{m}}}}=\left\|{{{\mathbf{w}}_{{D_{k}}}}-{{\mathbf{w}}_{{E_{m}}}}}\right\|$
signifies the Euclidean distance between $D_{k}$ and $E_{m}$,
$\xi_{{D_{k}}{E_{m}}}$ follows an exponential distribution with unit mean, and
$\alpha>2$ denotes the path-loss exponent. Similarly, the distance between
$D_{k}$ and $E_{m}$ is expressed as

$\displaystyle{d_{{D_{k}}{E_{m}}}}$
$\displaystyle=\left\|{{{\mathbf{w}}_{{D_{k}}}}-{{\mathbf{w}}_{{E_{m}}}}}\right\|$
(4)
$\displaystyle\geq\left|{\left\|{{{\mathbf{w}}_{{D_{k}}}}-{{{\mathbf{\hat{w}}}}_{{E_{m}}}}}\right\|-\left\|{{{{\mathbf{\hat{w}}}}_{{E_{m}}}}-{{\mathbf{w}}_{{E_{m}}}}}\right\|}\right|$
$\displaystyle\geq\left|{\underbrace{\left\|{{{\mathbf{w}}_{{D_{k}}}}-{{{\mathbf{\hat{w}}}}_{{E_{m}}}}}\right\|}_{\buildrel\Delta\over{=}{{\hat{d}}_{{D_{k}}{E_{m}}}}}-{r_{{E_{m}}}}}\right|.$

where ${{\hat{d}}_{{D_{k}}{E_{m}}}}$ signifies the estimation of
${d_{{D_{k}}{E_{m}}}}$. Then ${{h_{SE_{m}}}}$ and ${h_{{D_{k}}{E_{m}}}}$ are
approximated by

${{\hat{h}}_{S{E_{m}}}}\left(n\right)=\frac{{{{\beta_{0}}}}}{{{{\left({{{\hat{d}}_{S{E_{m}}}}\left(n\right)-{r_{{E_{m}}}}}\right)}^{2}}+{H^{2}}}},$
(5)

and

${{\hat{h}}_{{D_{k}}{E_{m}}}}=\frac{{{{\beta_{0}}}}}{{{{A_{k,m}}^{\alpha}}}}{\xi_{{D_{k}}{E_{m}}}},$
(6)

respectively, where
${A_{k,m}}=\left|{{{\hat{d}}_{{D_{k}}{E_{m}}}}-{r_{{E_{m}}}}}\right|$.

Due to the presence of FD eavesdropper nodes, the signal-to-interference-plus-
noise ratio (SINR) at $D_{k}$ is expressed as

${{\gamma_{{D_{k}}}}\left(n\right)=\frac{{{P_{S}}\left(n\right){h_{S{D_{k}}\left(n\right)}}}}{{\sum\limits_{i=1}^{M}{{P_{{E_{m}}}}\left(n\right){\hat{h}_{{D_{k}}{E_{m}}}}}+{\sigma^{2}}}},}$
(7)

where ${{P_{S}}\left(n\right)}$ and ${P_{{E_{m}}}}\left(n\right)$ denote the
transmit power of the ${S}$ and $E_{m}$, respectively, and ${{\sigma^{2}}}$
signifies the variance of additive white Gaussian noise (AWGN).

With the coordinated multi-point (CoMP) reception, all the cognitive users
jointly decode the received legitimate messages via maximal ratio combining
(MRC) scheme [29]. Then the achievable rate at the legitimate receiver is
expressed as

$\displaystyle{R_{S}}\left(n\right)$
$\displaystyle\approx{{\mathbb{E}}_{\xi}}\left[{{{\log}_{2}}\left({1+\sum\limits_{k=1}^{K}{{\gamma_{{D_{k}}}}}\left(n\right)}\right)}\right]$
(8)
$\displaystyle={{\mathbb{E}}_{\xi}}\left[{{{\log}_{2}}\left({1+\sum\limits_{k=1}^{K}{\frac{{{P_{S}}\left(n\right){h_{S{D_{k}}}}\left(n\right)}}{{\sum\limits_{m=1}^{M}{{P_{{E_{m}}}}\left(n\right){\hat{h}_{{D_{k}}{E_{m}}}}}+{\sigma^{2}}}}}}\right)}\right],$

where ${\mathbb{E}}\left[\cdot\right]$ denotes the expectation.

Since ${\hat{h}_{{D_{k}}{E_{m}}}}\left(n\right)$ in (6) is convex with respect
to $\xi$, ${R_{S}}\left(n\right)$ in (8) is convex with respect to
${{\hat{h}}_{{D_{k}}{E_{m}}}}\left(n\right)$, based on Jensen’s inequality, we
approximate the achievable at the legitimate receiver as [33]

$\displaystyle{R_{S}}\left(n\right)$
$\displaystyle\geq{\log_{2}}\left({1+\sum\limits_{k=1}^{K}{\frac{{{P_{S}}\left(n\right){h_{S{D_{k}}}\left(n\right)}}}{{\sum\limits_{m=1}^{M}{{P_{{E_{m}}}}\left(n\right){{\mathbb{E}}_{\xi}}\left[{{\hat{h}_{{D_{k}}{E_{m}}}}}\right]}+{\sigma^{2}}}}}}\right)$
(9)
$\displaystyle={\log_{2}}\left({1+\sum\limits_{k=1}^{K}{\frac{{{P_{S}}\left(n\right){h_{S{D_{k}}}\left(n\right)}}}{{{\beta_{0}}}{\sum\limits_{m=1}^{M}{{P_{{E_{m}}}}\left(n\right){{A_{k,m}^{-\alpha}}}}+{\sigma^{2}}}}}}\right)$
$\displaystyle\buildrel\Delta\over{=}R_{S}^{\mathbf{L}}\left(n\right),$

where the superscript ‘L’ signifies the lower bound.

Similarly to (7), the SINR at the $m$th $E$ can also be expressed as

${{\gamma_{S{E_{m}}}}\left(n\right)=\frac{{{P_{S}}\left(n\right){{\hat{h}}_{S{E_{m}}}\left(n\right)}}}{{\sum\limits_{j\neq
m}{{P_{{E_{j}}}}\left(n\right){h_{{E_{j}}{E_{m}}}}}+{\sigma^{2}}}},}$ (10)

where ${{h_{{E_{j}}{E_{m}}}}}$ denotes the interference from other
eavesdroppers.

It is assumed that all the eavesdroppers transmit signals with its maximum
power over the whole period and cooperatively intercept/decode the
confidential message from $S$. Thus, we have
${P_{{E_{m}}}}\left(n\right)={P_{E}}$. Moreover, the worst-case scenario is
considered, the interference from other eavesdroppers is assumed to be
perfectly cancelled [20]. Therefore, we obtain the SINR at the illegitimate
receiver as

$\displaystyle{R_{E}}\left(n\right)$
$\displaystyle={\log_{2}}\left({1+\sum\limits_{m=1}^{M}{{\gamma_{{SE_{m}}}}}}\right)$
(11)
$\displaystyle={\log_{2}}\left({1+\sum\limits_{m=1}^{M}{\frac{{{P_{S}}\left(n\right){{\hat{h}}_{S{E_{m}}}\left(n\right)}}}{{{P_{E}}\sum\limits_{j\neq
m}{{h_{{E_{j}}{E_{m}}}}}+{\sigma^{2}}}}}}\right)$
$\displaystyle={\log_{2}}\left({1+\sum\limits_{m=1}^{M}{\frac{{{P_{S}}\left(n\right){{\hat{h}}_{S{E_{m}}}\left(n\right)}}}{{{\sigma^{2}}}}}}\right).$

Then the achievable secrecy rate of the considered system is approximated as

$\displaystyle{R_{\sec}}\left(n\right)$
$\displaystyle={\left[{{{\log}_{2}}\left({R_{S}^{\mathbf{L}}\left(n\right)-R_{E}\left(n\right)}\right)}\right]^{+}}$
(12)
$\displaystyle={\left[{{{\log}_{2}}\left({\frac{{1+\sum\limits_{k=1}^{K}{\frac{{{P_{S}}\left(n\right){h_{S{D_{k}}}\left(n\right)}}}{{{{\beta_{0}}}{P_{E}}\sum\limits_{m=1}^{M}{{{A_{k,m}^{-\alpha}}}}+{\sigma^{2}}}}}}}{{1+\sum\limits_{m=1}^{M}{\frac{{{P_{S}}\left(n\right){{\hat{h}}_{S{E_{m}}}\left(n\right)}}}{{{\sigma^{2}}}}}}}}\right)}\right]^{+}},$

where ${\left[x\right]^{+}}=\max\left\\{{0,x}\right\\}$. Then the SEE of the
considered system is expressed as [23]

$\mu=\frac{{\sum\limits_{n=1}^{N}{{R_{\sec}}\left(n\right)}}}{{\sum\limits_{n=1}^{N}{{P_{S}}\left(n\right)}}}.$
(13)

At the $n$ time slot, the interference at $U_{r}$ is expressed as

$\displaystyle{{\hat{I}}_{r}}\left(n\right)$
$\displaystyle={P_{S}}\left(n\right){h_{S{U_{r}}}}\left(n\right)+{P_{E}}\sum\limits_{m=1}^{M}{{\mathbb{E}_{\xi}}\left[{{h_{{E_{m}}{U_{r}}}}}\right]}\hfill$
(14)
$\displaystyle\mathop{\leq}\limits^{\left(a\right)}{P_{S}}\left(n\right){h_{S{U_{r}}}}\left(n\right)+{P_{E}}\sum\limits_{m=1}^{M}{\frac{{{\beta_{0}}}}{{d_{{E_{m}}{U_{r}}}^{\alpha}}}}$
$\displaystyle\triangleq{I_{r}}\left(n\right)\hfill,$

where
${h_{{E_{m}}{U_{r}}}}=\frac{{{\beta_{0}}}}{{d_{{E_{m}}{U_{r}}}^{\alpha}}}{\xi_{{E_{m}}{U_{r}}}}$
denotes the channel coefficient between ${E_{m}}$ and ${U_{r}}$,
$\xi_{{E_{m}}{U_{r}}}$ follows an exponential distribution with unit mean,
step ${\left(a\right)}$ is obtained based on Jensen’s inequality [34] and
${d_{{E_{m}}{U_{r}}}}=\left|{\left\|{{{{\mathbf{w}}}_{{U_{r}}}}\left(n\right)-{{{\mathbf{\hat{w}}}}_{{E_{m}}}}}\right\|-{r_{{E_{m}}}}}\right|$.
To ensure the quality of service (QoS) of the primary users, the average
interference caused by $S$ and $E_{m}$ must be limited to the interference
temperature (IT) threshold, ${\Gamma_{r}}$, [19], [25], and [26]. Thus, we
have

$\frac{1}{N}\sum\limits_{n=1}^{N}{{I_{r}}\left(n\right)}\leq{\Gamma_{r}},\forall
r.$ (15)

##  III Problem Formulation

This work optimizes the WASR of the aerial CRNs by designing $S$’s horizontal
trajectory and transmission power, while the SEE is guaranteed. Let
${\mathbf{Q}}=\left\\{{{{\mathbf{q}}_{S}}\left(n\right),\forall n}\right\\}$,
${\mathbf{P}}=\left\\{{{P_{S}}\left(n\right),\forall n}\right\\}$. Thus, the
optimization problem is formulated as

$\displaystyle\mathcal{P}_{1}\,:\,$
$\displaystyle\mathop{\max}\limits_{{\mathbf{P}},{\mathbf{Q}}}\frac{1}{N}\sum\limits_{n=1}^{N}{R_{\sec}^{1}\left(n\right)}$
(16a) $\displaystyle{\mathrm{s.t.}}\;$ $\displaystyle{\mu\geq\Psi,}$ (16b)
$\displaystyle 0\leq{P_{S}}\left(n\right)\leq P_{S}^{\max},\forall n,$ (16c)
$\displaystyle\frac{1}{N}\sum\limits_{n=1}^{N}{{I_{r}}\left(n\right)}\leq{\Gamma_{r}},\forall
r,$ (16d)
$\displaystyle{{\mathbf{q}}_{S}}\left(0\right)={{\mathbf{q}}_{S}^{\mathbf{I}}},{{\mathbf{q}}_{S}}\left(n\right)={{\mathbf{q}}_{S}^{\mathbf{F}}},$
(16e)
$\displaystyle\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{q}}_{S}}\left({n-1}\right)}\right\|\leq
V_{S}^{\max}{\delta_{t}},\forall n,$ (16f)

where $\Psi$ is the SEE threshold, $P_{S}^{\max}$ is the peak transmit power
of $S$, ${\mathbf{q}}_{S}^{\mathbf{I}}$ and ${\mathbf{q}}_{S}^{\mathbf{F}}$
signify the initial and final position of $S$, respectively, $V_{S}^{\max}$
denotes the maximum speed of $S$, (16b) is the SEE constraint, (16c) is the
power constraint, (16d) is the constraint to ensure that the primary users
work properly, (16e) is the constraint for the initial and final position of
UAVs, and (16f) is the constraint for the maximum flight path in adjacent
periods of flight.

$\mathcal{P}_{1}$ is challengeable to solve since the operator
${\left[x\right]^{+}}$ makes the objective function non-smooth at zero value.
Moreover, $\mathcal{P}_{1}$ is a multivariate coupled optimization problem,
which is non-convex concerning horizontal trajectory variables ${\mathbf{Q}}$
and the transmit power variables ${\mathbf{P}}$.

##  IV Proposed Algorithm for Problem $\mathcal{P}_{1}$

To deal with the non-smoothness of ${R_{\sec}}\left(n\right)$, according to
Lemma 1 in [18] and [19], removing the operator ${\left[x\right]^{+}}$ does
not affect the optimal value of $\mathcal{P}_{1}$ because the secrecy rate can
be obtained by setting the transmit power to zero when it is less than zero.
Simultaneously, based on the alternating optimization (AO) method,
${{\mathbf{Q}}}$ and ${{\mathbf{P}}}$ are optimized in an alternating manner,
taking into account the other given variables.

###  IV-A Subproblem 1: Optimization Transmit Power of $S$

In this subsection, we optimize ${\mathbf{P}}$, the transmit power of $S$,
with given ${\mathbf{Q}}$. The original problem $\mathcal{P}_{1}$ is rewritten
as

$\displaystyle\mathcal{P}_{1.1}\,:\,$
$\displaystyle\max\limits_{{{\mathbf{P}}}}\;\frac{1}{N}\sum\limits_{n=1}^{N}{R_{\sec}^{1}\left(n\right)}$
(17a) $\displaystyle{\mathrm{s.t.}}\;$ $\displaystyle\;{\mu_{1}\geq\Psi,}$
(17b) $\displaystyle\;(\textrm{\ref{P1c}}),(\textrm{\ref{P1d}}),$ (17c)

where
${\mu_{1}}=\frac{{\sum\limits_{n=1}^{N}{R_{\sec}^{1}\left(n\right)}}}{{\sum\limits_{n=1}^{N}{{P_{S}}\left(n\right)}}}$,
$R_{\sec}^{1}\left(n\right)={{\log}_{2}}\left({1+{A_{n}}{P_{S}}\left(n\right)}\right)-{{\log}_{2}}\left({1+{B_{n}}{P_{S}}\left(n\right)}\right)$,
${A_{n}}=\sum\limits_{k=1}^{K}{\frac{{{h_{S{D_{k}}}\left(n\right)}}}{{{{\beta_{0}}}{P_{E}}\sum\limits_{m=1}^{M}{{{A_{k,m}^{-\alpha}}}}+{\sigma^{2}}}}}$,
and
${B_{n}}=\sum\limits_{i=1}^{M}{\frac{{{\hat{h}_{S{E_{m}}}}\left(n\right)}}{{{\sigma^{2}}}}}$.

To tackle the non-convex characterization in (17a) and (17b), the first-order
Taylor expansion is utilized to approximate $R_{\sec}^{1}\left(n\right)$ as

$\displaystyle R_{\sec}^{2}\left(n\right)$
$\displaystyle={{\log}_{2}}\left({1+{A_{n}}{P_{S}}\left(n\right)}\right)-{{\log}_{2}}\left({1+{B_{n}}P_{S}^{\left(l\right)}\left(n\right)}\right)$
(18)
$\displaystyle-\frac{{{B_{n}}}}{{\ln\left(2\right)\left({1+{B_{n}}P_{S}^{\left(l\right)}\left(n\right)}\right)}}\left({{P_{S}}\left(n\right)-P_{S}^{\left(l\right)}\left(n\right)}\right),$

where $P_{S}^{\left(l\right)}\left(n\right)$ is a given feasible point in the
${l}$th iteration. Then $\mathcal{P}_{1.1}$ is approximated as

$\displaystyle\mathcal{P}_{1.2}\,:\,$
$\displaystyle\max\limits_{{{\mathbf{P}}}}\;\frac{1}{{N}}\sum\limits_{n=1}^{N}{R_{\sec}^{2}\left(n\right)}$
(19a) $\displaystyle{\mathrm{s.t.}}\;$
$\displaystyle\;\sum\limits_{n=1}^{N}{R_{\sec}^{2}\left(n\right)}\geq\Psi\sum\limits_{n=1}^{N}{{P_{S}}\left(n\right)},$
(19b) $\displaystyle\;(\textrm{\ref{P1c}}),(\textrm{\ref{P1d}}).$ (19c)

$\mathcal{P}_{1.2}$ is a convex problem and can be solved by existing
optimization tools such as CVX.

###  IV-B Subproblem 2: Optimizing Trajectory of $S$

Given ${\mathbf{P}}$, $\mathcal{P}_{1}$ is rewritten as

$\displaystyle\mathcal{P}_{2.1}\,:\,$
$\displaystyle\;\max\limits_{{{\mathbf{Q}}}}\;{\frac{1}{N}\sum\limits_{n=1}^{N}{R_{\sec}^{3}\left(n\right)}}$
(20a) $\displaystyle{\mathrm{s.t.}}\;$ $\displaystyle\;{\mu_{2}\geq\Psi,}$
(20b)
$\displaystyle\;\frac{1}{N}\sum\limits_{n=1}^{N}{I_{r}^{\mathbf{Q}}\left(n\right)}\leq{\Gamma_{r}},\forall
r,$ (20c) $\displaystyle\;(\textrm{\ref{P1e}}),(\textrm{\ref{P1f}}),$ (20d)

where
${\mu_{2}}=\frac{{\sum\limits_{n=1}^{N}{R_{\sec}^{3}\left(n\right)}}}{{\sum\limits_{n=1}^{N}{{P_{S}}\left(n\right)}}}$,
$R_{\sec}^{3}\left(n\right)={\log_{2}}\left({1+\sum\limits_{k=1}^{K}{\frac{{{f_{1,k}}}}{{d_{S{D_{k}}}^{2}+{H^{2}}}}}}\right)$
$-{\log_{2}}\left({1+\sum\limits_{m=1}^{M}{\frac{{{f_{2}}}}{{{{\left({{{\hat{d}}_{S{E_{m}}}}-{r_{{E_{m}}}}}\right)}^{2}}+{H^{2}}}}}}\right)$,
${f_{1,k}}=\frac{{{P_{S}}\left(n\right)}}{{{P_{E}}\sum\limits_{m=1}^{M}{{A_{k,m}^{-\alpha}}}+\frac{1}{{{\rho_{0}}}}}}$,
${f_{2}}={\rho_{0}}{P_{S}}\left(n\right)$,
$I_{r}^{\mathbf{Q}}\left(n\right)=\frac{{{P_{S}}\left(n\right){\rho_{0}}}}{{{d_{S{U_{r}}}^{2}}+{H^{2}}}}+{P_{E}}\sum\limits_{m=1}^{M}{h_{{E_{m}}{U_{r}}}}$,
and ${\rho_{0}}=\frac{{{\beta_{0}}}}{{{\sigma^{2}}}}$.

It must be noted that $\mathcal{P}_{2.1}$ is a non-convex problem because the
objective function in (20a) is a non-concave function with respect to
${{{\mathbf{q}}_{S}}\left(n\right)}$ and (20b) and (20c) are non-convex
constraints with respect to ${{{\mathbf{q}}_{S}}\left(n\right)}$.

By introducing new slack variables and using the successive convex
approximation (SCA), we obtain

$\displaystyle{{R_{\sec}^{3}\left(n\right)}}$
$\displaystyle\geq{\log_{2}}\left({1+\sum\limits_{k=1}^{K}{\frac{{{f_{1,k}}}}{{{\tilde{\zeta}_{k}}^{\left(l\right)}\left(n\right)}}}}\right)$
(21)
$\displaystyle-{\log_{2}}\left({1+\sum\limits_{m=1}^{M}{\frac{{{f_{2}}}}{{{{\tilde{\xi}_{m}}}\left(n\right)}}}}\right)$
$\displaystyle-\frac{{\sum\limits_{k=1}^{K}{\left({\frac{{{f_{1,k}}}}{{{{\left({{\tilde{\zeta}_{k}}^{\left(l\right)}\left(n\right)}\right)}^{2}}}}\left({{{\tilde{\zeta}_{k}}}\left(n\right)-{\tilde{\zeta}_{k}}^{\left(l\right)}\left(n\right)}\right)}\right)}}}{{\ln\left(2\right)\left({1+\sum\limits_{k=1}^{K}{\frac{{{f_{1,k}}}}{{{\tilde{\zeta}_{k}}^{\left(l\right)}\left(n\right)}}}}\right)}}$
$\displaystyle\buildrel\Delta\over{=}R_{\sec}^{3,{\mathbf{L}}}\left(n\right),$

where the superscript ‘$l$’ denotes the iteration index,
${{\tilde{\zeta}_{k}}}\left(n\right)$ and ${{\tilde{\xi}_{m}}}\left(n\right)$
are new slack variables, which must satisfy

${{\tilde{\zeta}_{k}}}\left(n\right)\geq{d_{S{D_{k}}}^{2}}+{H^{2}}$ (22)

and

${{\tilde{\xi}_{m}}}\left(n\right)\leq{\left({{{\hat{d}}_{S{E_{m}}}}-{r_{{E_{m}}}}}\right)^{2}}+{H^{2}},$
(23)

respectively. Due to the right-hand-sight of (23) is convex which makes (23)
become non-convex. According to SCA, (23) is rewritten as

$\displaystyle{{\tilde{\xi}_{m}}}\left(n\right)$
$\displaystyle\leq{H^{2}}-2{r_{{E_{m}}}}\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{{E_{m}}}}}\right\|+r_{{E_{m}}}^{2}$
(24)
$\displaystyle+2{\left({{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{{E_{m}}}}}\right)^{T}}\left({{{\mathbf{q}}_{S}}\left(n\right)-{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)}\right)$
$\displaystyle+{\left\|{{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)-{{\mathbf{w}}_{{E_{m}}}}}\right\|^{2}}.$

With the same method, (20c) is rewritten as

$\frac{1}{N}\sum\limits_{n=1}^{N}{\left({\left.{\frac{{{\rho_{0}}{P_{S}}\left(n\right)}}{{{{\tilde{d}}_{S{U_{r}}}}+{H^{2}}}}+{P_{E}}\sum\limits_{m=1}^{M}{h_{{E_{m}}{U_{r}}}}}\right)}\right.}\leq{\Gamma_{r}},$
(25)

where ${{{\tilde{d}}_{S{U_{r}}}}}$ is a new slack variable which satisfies

$\displaystyle{{\tilde{d}}_{S{U_{r}}}}$
$\displaystyle\leq{\left\|{{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)-{{\mathbf{w}}_{r}}}\right\|^{2}}$
(26)
$\displaystyle+2{\left({{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{r}}}\right)^{T}}\left({{{\mathbf{q}}_{S}}\left(n\right)-{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)}\right).$

Based on the above conversion, $\mathcal{P}_{2.1}$ is approximated as

$\displaystyle{\mathcal{P}_{2.2}}$
$\displaystyle\mathop{\max}\limits_{{{\mathbf{q}}_{S}},{{\tilde{\xi}_{m}}}\left(n\right),{{\tilde{\zeta}_{k}}}\left(n\right),{{\tilde{d}}_{S{U_{r}}}}}{{\frac{1}{N}\sum\limits_{n=1}^{N}{R_{\sec}^{3,{\mathbf{L}}}\left(n\right)}}}$
(27a) $\displaystyle{\mathrm{s.t.}}\;$ $\displaystyle\;{\mu_{3}\geq\Psi,}$
(27b)
$\displaystyle\frac{1}{N}\sum\limits_{n=1}^{N}{\frac{{{\beta_{0}}{P_{S}}\left(n\right)}}{{{{\tilde{d}}_{S{U_{r}}}}+{H^{2}}}}}+{P_{E}}\sum\limits_{m=1}^{M}{h_{{E_{m}}{U_{r}}}}\leq{\Gamma_{r}},\forall
r,$ (27c)
$\displaystyle{{\tilde{\zeta}_{k}}}\left(n\right)\geq{d_{S{D_{k}}}^{2}}+{H^{2}},\forall
n,\forall k,$ (27d)
$\displaystyle{{\tilde{\xi}_{m}}}\left(n\right)\leq{H^{2}}-2{r_{{E_{m}}}}\left\|{{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{{E_{m}}}}}\right\|+r_{{E_{m}}}^{2}$
$\displaystyle\;\;\;\;\;\;\;\;+2{\left({{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{{E_{m}}}}}\right)^{T}}\left({{{\mathbf{q}}_{S}}\left(n\right)-{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)}\right)$
(27e)
$\displaystyle\;\;\;\;\;\;\;\;+{\left\|{{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)-{{\mathbf{w}}_{{E_{m}}}}}\right\|^{2}},\forall
n,\forall m,$
$\displaystyle{{\tilde{d}}_{S{U_{r}}}}\leq{\left\|{{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)-{{\mathbf{w}}_{r}}}\right\|^{2}}$
$\displaystyle\;\;+2{\left({{{\mathbf{q}}_{S}}\left(n\right)-{{\mathbf{w}}_{r}}}\right)^{T}}\left({{{\mathbf{q}}_{S}}\left(n\right)-{\mathbf{q}}_{S}^{\left(l\right)}\left(n\right)}\right),\forall
n,\forall r,$ (27f) $\displaystyle(\textrm{\ref{P1e}}),(\textrm{\ref{P1f}}),$
(27g)

where
${\mu_{3}}=\frac{{\sum\limits_{n=1}^{N}{R_{\sec}^{3,{\mathbf{L}}}\left(n\right)}}}{{\sum\limits_{n=1}^{N}{{P_{S}}\left(n\right)}}}$.
Now $\mathcal{P}_{2.2}$ is convex and can be solved with existing optimisation
tools.

TABLE III: List of Simulation Parameters. Parameter | Value
---|---
${\mathbf{q}}_{S}^{0}$ |  ${\left[{-300,-200}\right]^{T}}$, ${\left[{-250,800}\right]^{T}}$
${\mathbf{q}}_{S}^{F}$ |  ${\left[{100,-200}\right]^{T}}$, ${\left[{0,-200}\right]^{T}}$
${{\mathbf{w}}_{{D_{k}}}}$ |  ${\left[{-200,500;-100,500;0,500}\right]^{T}}$ ${\left[{-100,300;-50,400;0,300}\right]^{T}}$ ${\left[{-250,300;-100,0;50,300}\right]^{T}}$
${{\mathbf{\hat{w}}}_{{E_{m}}}}$ |  ${\left[{-150,700;-50,700}\right]^{T}}$ ${\left[{-125,100;-150,200}\right]^{T}}$ ${\left[{-200,-20;0,-20}\right]^{T}}$
${{\mathbf{w}}_{{U_{r}}}}$ |  ${\left[{-100,0}\right]^{T}}$ ${\left[{-200,400}\right]^{T}}$ ${\left[{-100,400}\right]^{T}}$
${r_{{E_{m}}}}$ | 10 m, 30 m
$H$ | 100 m
${\Gamma_{r}}$ | -110 dBm
${\sigma^{2}}$ | -110 dBm
$P_{S}^{\max}$ | 3 W
${P_{{E}}}$ | 0.1 W
$\alpha$ | 2.2
$V_{S}^{\max}$ | 60 m/s
${\beta_{0}}$ | \- 60 dB
${\delta_{t}}$ | 1 s
$\varepsilon$ | 0.01
$\Psi$ | 1

###  IV-C Convergence Analysis of Algorithm 1

Combining the above two subproblems and applying BCD, an iterative algorithm
for solving $\mathcal{P}_{1}$ is proposed. The suboptimal solution is obtained
by selecting the initial feasible points and solving $\mathcal{P}_{1.2}$ and
$\mathcal{P}_{2.2}$ alternatively. Each iteration’s solutions serve as
feasible input values at the next iteration. The objective function of the
original problem $\mathcal{P}_{1}$ at the $l$th iteration is defined as
${R_{\sec}}\left({{{\mathbf{P}}^{l}},{\mathbf{Q}}^{l}}\right)$, meanwhile, the
details of the entire iteration of $\mathcal{P}_{1}$ are summarized in
Algorithm 1, where the tolerance of the convergence of the algorithm is
denoted as $\varepsilon$. The convergence of Algorithm 1 is proved as follows.

Input: Initialization of feasible points.

while
_${R_{\sec}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{\mathbf{Q}}^{\left(l+1\right)}}\right)-{R_{\sec}}\left({{{\mathbf{P}}^{\left(l\right)}},{\mathbf{Q}}^{\left(l\right)}}\right)\geq{\varepsilon}$_
do

1\. Solve ($\mathcal{P}_{1.2}$) for given ${{\mathbf{Q}}^{\left(l\right)}}$
and obtain the solution ${{\mathbf{P}}^{{\left(l+1\right)}}}$;

2\. Solve ($\mathcal{P}_{2.2}$) for given
${{{\mathbf{P}}^{{\left(l+1\right)}}}}$ and obtain the solution
${{\mathbf{Q}}^{{\left(l+1\right)}}}$;

3\. $l=l+1$;

4\. Calculate
${R_{\sec}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l+1\right)}}}\right)$.

end while

Output:
${R_{\sec}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l+1\right)}}}\right)$
with
${{\mathbf{P}}^{*}}={{\mathbf{P}}^{\left(l+1\right)}},\;{{\mathbf{Q}}^{*}}={{\mathbf{Q}}^{\left(l+1\right)}}$.

Algorithm 1 Iterative Procedure of Problem $\mathcal{P}_{1}$

###### Proof.

In Step 1 of Algorithm 1, the power allocation of UAV solved by applying
interior point method. Thus, we have

$\displaystyle{R_{\sec}}\left({{{\mathbf{P}}^{\left(l\right)}},{{\mathbf{Q}}^{\left(l\right)}}}\right)$
$\displaystyle\leq{{R_{\sec}^{\mathbf{L},1}}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l\right)}}}\right)$
(28)
$\displaystyle\leq{R_{\sec}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l\right)}}}\right),$

where the objective function of problem $\mathcal{P}_{1.2}$ is lower-bounded
by
${{R_{\sec}^{\mathbf{L},1}}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l\right)}}}\right)$.
Then, by solving $\mathcal{P}_{2.2}$ the suboptimal solution
${{{\mathbf{Q}}^{\left(l+1\right)}}}$ is obtained, and have

$\displaystyle{R_{\sec}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{l}}}\right)$
$\displaystyle\leq{{R_{\sec}^{\mathbf{L},2}}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l+1\right)}}}\right)$
(29)
$\displaystyle\leq{R_{\sec}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l+1\right)}}}\right),$

where the objective function of problem $\mathcal{P}_{2.2}$ is lower-bounded
by
${{R_{\sec}^{\mathbf{L},2}}}\left({{{\mathbf{P}}^{\left(l+1\right)}},{{\mathbf{Q}}^{\left(l+1\right)}}}\right)$.

The objective function is always non-decreasing after each iteration, as shown
in (29). Meanwhile, we additionally notice that due to the feasible set under
the constraints, the objective function in $\mathcal{P}_{1}$ is upper bounded
by a finite value. Therefore, Algorithm 1 is convergent. ∎

The complexity of Algorithm 1 comes from two aspects. One is optimizing
transmit power of $S$, the other is solving $S$’s trajectory. They are solved
by the interior point method with complexity ${\cal
O}\left({{\phi^{3.5}}\ln\left({\frac{1}{\varsigma}}\right)}\right)$ and ${\cal
O}\left({{\varphi^{3.5}}\ln\left({\frac{1}{\varsigma}}\right)}\right)$,
respectively, where $\phi=2N+R+1$, ${\varphi=NK+NM+NR+N+R+5}$, and $\varsigma$
denotes convergence precision. It is assumed that the number of iterations is
$I$. Thus, the total complexity of Algorithm 1 is ${\cal
O}\left({I\ln\left({\frac{1}{\varsigma}}\right)\left({{\phi^{3.5}}+{\varphi^{3.5}}}\right)}\right)$.

##  V Numerical Results and Analysis

(a) The optimal trajectory of $S$.

(b) The transmit power of $S$.

(c) The achievable secrecy rate.

Figure 2: Scenario 1 wherein the FD-eavesdroppers are uniformly distributed
with $K=3$ terrestrial users.

(a) The optimal trajectory of $S$.

(b) The transmit power of $S$.

(c) The achievable secrecy rate.

Figure 3: Scenario 2 wherein the FD-eavesdroppers are cluster distributed with
$K=3$ terrestrial users.

(a) The optimal trajectory of $S$.

(b) The transmit power of $S$.

(c) The achievable secrecy rate.

Figure 4: Scenario 3 wherein the FD-eavesdroppers are interweaved in space
with $K=3$ terrestrial users.

In this section, simulation results are presented to verify the performance of
the joint transmit power control and UAV trajectory optimization algorithm. To
verify the effectiveness of our proposed algorithm, three scenarios are
considered, in which the two FD-eavesdroppers are uniformly distributed,
cluster distributed, and interweaved distributed with three terrestrial users,
respectively. The details of parameter setup are listed in Table III. Similar
to [18], [19], [20], the following benchmark schemes are also considered : 1)
Benchmark 1: $S$ operates with a fixed transmission power and a fixed
trajectory, which is the initial value for the optimized transmit power and
trajectory (denoted as FPFT). 2) Benchmark 2: $S$ works with fixed transmit
power and the trajectory is optimized (denoted as FPOT). 3) Benchmark 3: $S$
works with fixed trajectory and its transmission power is optimized (denoted
as OPFT).

To illustrate the effect of colluding FD eavesdroppers, the trajectory
obtained with different schemes for the different scenarios is shown in Figs.
2 \- 4. For the FPOT scheme where $S$ works with fixed power, it should try to
stay away from the area where the eavesdroppers are located while getting
close to the legitimate users. For the OPFT scheme where $S$ flies with a
given trajectory, the transmit power should be reduced (to zero) as it
approaches the area where the eavesdroppers are located to avoid the
information being eavesdropped. In the proposed scheme, the power is gradually
increased/decreased to maximize the secrecy rate as it approaches the
legitimate users (such as $D_{1}$ in Figs. 2(a) and 4(a)) and ensure that the
interference to the primary user does not exceed the IT threshold. As $S$
flies, it gradually moves away from the legitimate users and approaches the
area where the eavesdroppers are located, and the transmission power is
reduced to avoid the information being wiretapped. The characteristics of FPOT
and OPFT are considered simultaneously to obtain the optimal secrecy
performance. It can be observed from Fig. 4(c) that the achievable secrecy
rate of FPOT is higher than that of the proposed at 31s. This is because the
optimized transmit power of the proposed is lower to decrease the interference
to the primary user. We note that the optimality of the proposed scheme does
not mean the proposed scheme is time-piecewise optimal but time-averaged
optimal.

Figure 5: The WASR versus the number of iterations for different scenarios.

Fig. 5 shows the convergence of the proposed scheme in terms of the WASR
versus the number of iterations in different scenarios. The results illustrate
the convergence of the proposed scheme in terms of jointly optimized transmit
power and trajectory for the aerial base station. One can observe that the
WASR increases quickly within a small number of iterations and converges
within around eight iterations in the listed three scenarios.

Figure 6: The WASR versus ${\Gamma_{r}}$ with different schemes.

Fig. 6 shows the effect of the IT threshold on the WASR of different schemes.
It can be observed that, as the threshold increases, the curves corresponding
to the proposed and OPFT schemes show an increasing trend while those to FPFT
and FPOT schemes remain constant. The reason is as ${\Gamma_{r}}$ increases,
the optimized transmit power increases, allowing the secrecy rate of the
system to be improved. In the lower-${\Gamma_{r}}$ region, the transmit power
of $S$ is lower since the average interference power caused by $S$ and
eavesdroppers must be limited to the IT threshold to ensure that the QoS of
the primary users is not affected. Thus, the secrecy rates with the proposed
scheme and FPOT scheme outperform those with FPFT and OPFT schemes. Moreover,
the curves for the proposed and OPFT schemes have the apparent increasing
trend because ${\Gamma_{r}}$ directly affects the transmit power of $S$ since
the transmit powers in the proposed and OPFT schemes are adaptively changed
while the transmit power in FPFT and FPOT schemes are fixed. When the IT
threshold reaches a certain value, the system becomes a non-cognitive mode
where the transmit power is limited by the maximum transmit power of $S$. Then
the WASR is independent of the IT threshold in the larger-${\Gamma_{r}}$
region.

Figure 7: The WASR of scenario 1 under different schemes versus $T$.

Fig. 7 indicates the WASR versus the flight period $T$ for the considered
system under different schemes. It can be observed that the WASR of the
considered system obtained by the proposed scheme outperforms those of others
schemes. This is because both the UAV’s trajectory and transmission power are
jointly optimized in the proposed scheme simultaneously.

##  VI Conclusion

This work investigated the security of aerial CRNs with multiple uncertainty
FD eavesdroppers. By jointly designing the aerial base station’s trajectory
and transmission power, the WASR is formulated as a non-convex maximization
problem. Based on BCD and SCA technologies, a new efficient iterative
algorithm is proposed to solve the non-convex problem, and suboptimal
solutions are obtained. Numerical results verified the convergence and
effectiveness of our proposed algorithm. The A2G links in this work are
assumed to be LoS links. In some scenarios, the A2G links are LoS links with
elevation angle-dependent and environment-dependent probability. Thus,
designing the 3D trajectory will be an interesting problem in the future work.
Considering the distribution of the eavesdroppers’ location is also an
interesting work and will be part of the future work.

## References

  * [1] Y. Zeng, Q. Wu, and R. Zhang, “Accessing from the sky: A tutorial on UAV communications for 5G and beyond,” Proc. IEEE, vol. 107, no. 12, pp. 2327-2375, Dec. 2019.
  * [2] B. Li, Z. Fei, and Y. Zhang, “UAV communications for 5G and beyond: Recent advances and future trends,” IEEE Internet Things J., vol. 6, no. 2, pp. 1-23, Apr. 2019.
  * [3] Q. Wu, L. Liu, and R. Zhang, “Fundamental trade-offs in communication and trajectory design for UAV-enabled wireless network,” IEEE Wireless Commun., vol. 26, no. 1, pp. 36-44, Feb. 2019.
  * [4] C. Yan, L. Fu, J. Zhang, and J. Wang, “A comprehensive survey on UAV communication channel modeling,” IEEE Access, vol. 7, pp. 107769-107792, Aug. 2019.
  * [5] Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for multi-UAV enabled wireless networks,” IEEE Trans. Wireless Commun., vol. 17, no. 3, pp. 2109-2121, Mar. 2018.
  * [6] A. Bejaoui, K.-H. Park, and M.-S. Alouini, “A QoS-oriented trajectory optimization in swarming unmanned-aerial-vehicles communications,” IEEE Wireless Commun. Lett., vol. 9, no. 6, pp. 791-794, Jun. 2020.
  * [7] J. Wang, C. Jiang, Z. Wei, C. Pan, H. Zhang, and Y. Ren, “Joint UAV hovering altitude and power control for space-air-ground IoT networks,” IEEE Internet Things J., vol. 6, no. 2, pp. 1741-1753, Apr. 2019.
  * [8] L. Yang, H. Yao, J. Wang, C. Jiang, A. Benslimane, and Y. Liu, “Multi-UAV-enabled load-balance mobile-edge computing for IoT networks,” IEEE Internet Things J., vol. 7, no. 8, pp. 6898-6908, Aug. 2020.
  * [9] T. Liu, M. Cui, G. Zhang, Q. Wu, X. Chu, and J. Zhang, “3D trajectory and transmit power optimization for UAV-enabled multi-link relaying systems,” IEEE Trans. Green Commun. Netw., vol. 5, no. 1, pp. 392-405, Mar. 2021.
  * [10] J.-H. Lee, K.-H. Park, Y.-C. Ko, and M.-S. Alouini, “Throughput maximization of mixed FSO/RF UAV-aided mobile relaying with a buffer,” IEEE Trans. Wireless Commun., vol. 20, no. 1, pp. 683-694, Jan. 2021.
  * [11] C. Wang, X. Chen, J. An, Z. Xiong, C. Xing, N. Zhao, and D. Niyato, “Covert communication assisted by UAV-IRS,” IEEE Trans. Commun., vol. 71, no. 1, pp. 357-369, 2022.
  * [12] C. Zhan, Y. Zeng, and R. Zhang, “Energy-efficient data collection in UAV enabled wireless sensor network,” IEEE Wireless Commun. Lett., vol. 7, no. 3, pp. 328-331, Jun. 2018.
  * [13] W. Lu, P. Si, F. Lu, B. Li, Z. Liu, S. Hu, and Y. Gong, “Resource and trajectory optimization in UAV-powered wireless communication system,” Sci. China Inf. Sci., vol. 64, no. 4, pp. 1-14, Mar. 2021.
  * [14] W. Luo, Y. Shen, B. Yang, S. Wang, and X. Guan, “Joint 3-D trajectory and resource optimization in multi-UAV-enabled IoT networks wth wireless power transfer,” IEEE Internet Things J., vol. 8, no. 10, pp. 7833-7848, May. 2021.
  * [15] Q. Wu, W. Mei, and R. Zhang, “Safeguarding wireless network with UAVs: A physical layer security perspective,” IEEE Wireless Commun., vol. 26, no. 5, pp. 12-18, Oct. 2019.
  * [16] H.-M. Wang, X. Zhang, and J.-C. Jiang, “UAV-involved wireless physical-layer secure communications: Overview and research directions,” IEEE Wireless Commun., vol. 26, no. 5, pp. 32-39, Oct. 2019.
  * [17] Z. Ullah, F. Al-Turjman, and L. Mostarda, “Cognition in UAV-aided 5G and beyond communications: A survey,” IEEE Trans. Cogn. Commun. Netw., vol. 6, no. 3, pp. 872-891, Sept 2020.
  * [18] G. Zhang, Q. Wu, M. Cui, and R. Zhang, “Securing UAV communications via joint trajectory and power control,” IEEE Trans. Wireless Commun., vol. 18, no. 2, pp. 1376-1389, Feb. 2019.
  * [19] Y. Wang, L. Chen, Y. Zhou, X. Liu, F. Zhou, and N. Al-Dhahir, “Resource allocation and trajectory design in UAV-assisted jamming wideband cognitive radio networks,” IEEE Trans. Cogn. Commun. Netw., vol. 7, no. 2, pp. 635-647, Jun. 2021.
  * [20] B. Duo, J. Luo, Y. Li, H. Hu, and Z. Wang, “Joint trajectory and power optimization for securing UAV communications against active eavesdropping,” China Commun., vol. 18, no. 1, pp. 88-99, Jan. 2021.
  * [21] M. Cui, G. Zhang, Q. Wu, and D. W. K. Ng, “Robust trajectory and transmit power design for secure UAV communications,” IEEE Trans. Veh. Technol.,vol. 67, no. 9, pp. 9042-9046, Sept. 2018.
  * [22] C. Zhong, J. Yao, and J. Xu, “Secure UAV communication with cooperative jamming and trajectory control,” IEEE Commun. Lett., vol. 23, no. 2, pp. 286-289, Feb. 2019.
  * [23] R. Zhang, X. Pang, W. Lu, N. Zhao, Y. Chen, and D. Niyato, “Dual-UAV enabled secure data collection with propulsion limitation,” IEEE Trans. Wireless Commun., vol. 20, no. 11, pp. 7445-7459, Jun. 2021.
  * [24] W. Wang, X. Li, R. Wang, K. Cumanan, W. Feng, Z. Ding, and O. A. Dobre, “Robust 3D-Trajectory and time switching optimization for dual-UAV-enabled secure communications,” IEEE J. Sel. Areas Commun., vol. 39, no. 11, pp. 3334-3347, Nov. 2021.
  * [25] Y. Zhou, F. Zhou, H. Zhou, D. W. K. Ng, and R. Q. Hu, “Robust trajectory and transmit power optimization for secure UAV-enabled cognitive radio networks,” IEEE Trans. Commun., vol. 68, no. 7, pp. 4022-4034, Jul. 2020.
  * [26] P. X. Nguyen, V.-D. Nguyen, H. V. Nguyen, and O.-S. Shin, “UAV-assisted secure communications in terrestrial cognitive radio networks: Joint power control and 3D trajectory optimization,” IEEE Trans. Veh. Technol., vol. 70, no. 4, pp. 3298-3313, Apr. 2021.
  * [27] S. Liu, Y. Yu, L. Guo, P. L. Yeoh, B. Vucetic, Y. Li, and T. Q. Duong, “Satisfaction-maximized secure computation offloading in multi-eavesdropper MEC networks,” IEEE Trans. Wireless Commun., vol. 21, no. 6, pp. 4227-4241, Jun. 2022.
  * [28] H. Lei, H. Yang, K. H. Park, I. S. Ansari, J. Jiang, and M. S. Alouini, “Joint trajectory design and user scheduling for secure aerial underlay IoT systems,” IEEE Internet Things J.,vol. 10, no. 15, pp. 13637-13648, Aug. 2023.
  * [29] J. Yao and J. Xu, “Joint 3D maneuver and power adaptation for secure UAV communication With CoMP reception,” IEEE Trans. Wireless Commun., vol. 19, no. 10, pp. 6992-7006, Oct. 2020.
  * [30] Y. Gao, H. Tang, B. Li, and X. Yuan, “Energy minimization for robust secure transmission in UAV networks with multiple colluding eavesdroppers,” IEEE Commun. Lett., vol. 25, no. 7, pp. 2353-2357, Jul. 2021.
  * [31] H. Fu, Z. Sheng, A. A. Nasir, A. H. Muqaibel, and L. Hanzo, “Securing the UAV-aided non-orthogonal downlink in the face of colluding eavesdroppers,” IEEE Trans. Veh. Technol., vol. 71, no. 6, pp. 6837-6842, Jun. 2022.
  * [32] B. Liu, Y. Wan, F. Zhou, Q. Wu, and R. Q. Hu, “Resource allocation and trajectory design for MISO UAV-assisted MEC networks,” IEEE Trans. Veh. Technol., vol. 71, no. 5, pp. 4933-4948, May 2022.
  * [33] M. Hua, L. Yang, Q. Wu, and A. L. Swindlehurst, “3D UAV trajectory and communication design for simultaneous uplink and downlink transmission,” IEEE Trans. Veh. Technol., vol. 68, no. 9, pp. 5908-5923, Sept. 2020.
  * [34] S. Boyd and L. Vandenberghe, Convex Optimization: Cambridge Univ Press, 2004.