# Multipath Stability Routing in Cognitive UAV Swarm for Emergency Communications: A Hypergraph Matching Approach

Authors: Xiaozheng Ma; Yao Wang; Long Zhang
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

With the rapid development of unmanned aerial vehicles (UAVs), it has been considered as an effective solution for emergency communications. In order to solve the contradiction between the rapid growth of user equipments and the shortage of spectrum resources as well as the problem that a single UAV cannot meet the multimission requirements, we integrate the cognitive radio spectrum access and UAVs to form a cognitive UAV swarm. Since most of the base stations cannot work properly in disaster scenarios, cognitive UAV swarm transmits the collected information through multihop routing in the form of a decentralized distributed ad hoc network. In order to solve the problems of low robustness of single paths and unstable links, we explore a multipath stability optimization problem and then propose a multipath stable routing algorithm based on the multidimensional hypergraph matching method. The results show that compared with the greedy algorithm, the proposed algorithm can obtain higher path stability while obtaining multiple node disjoint paths. Finally, this paper provides a reference value for establishing multipath routing of UAV ad hoc networks in emergency communication scenarios.


---

## Body

## 1. Introduction

Various natural disasters occur from time to time around the world, which have caused huge economic losses and casualties. Disaster environments often disrupt roads and communications infrastructure. And due to the complex environment, things are unpredictable, and every minute of delay can lead to irreversible damage. People in disaster areas need to communicate with the outside world for assistance. The outside world also urgently needs to know the real-time situation of the disaster area in order to carry out rescue. Therefore, it is very necessary for the establishment of post-disaster emergency communications for communications and rescue in the disaster area.Limited by space and environment, the existing traditional methods lack flexibility and typically require central infrastructure such as base stations to provide support. To address these challenges, UAVs have recently attracted much attention for assisting communications due to their flexibility [1–3]. And the fact shows that the centralized infrastructure is not required through the form of the flying UAV ad hoc network and the cooperation between the swarm makes the processing more efficient [4, 5].Moreover, under the ad hoc network, using the traditional static access spectrum is no longer applicable, and the spectrum resources of other frequency bands are already crowded. Fortunately, we can dynamically access licensed subchannels on resource block through underlay spectrum access by building a cognitive radio module on the drone.Affected by the dynamics and instability of UAVs, single-path routing is easily interrupted when the sensing data is transmitted between UAVs, causing additional latency on retransmissions. To facilitate the distribution of the traffic load and also to avoid the route failure, there is a need to exploit the multipath routing from the perspective of path diversity to improve the multihop transmission efficiency and stability in UAV swarm ad hoc networks [6]. From the appropriate multipath routing perspective, it is expected that the stable multipath routing will increase both the end-to-end throughput and the connection reliability.We need to consider the following key points when designing the system: (1) how to design routing metric considering UAV dynamics and instability in a distributed manner without global information and (2) how to build node disjoint multipath stable paths.Based on the above requirements, we can find the importance of multipath routing combined with cognitive radio, the complexity of the environment, and the dynamics of the UAV to design stable routing. In this paper, we use hypergraph matching theory to find node-disjoint multipath routing. And according to our research, the problem of finding multipath routing through hypergraph matching has not been thoroughly researched in the relevant literature. To fill this gap, we establish a multipath stable routing based on hypergraph matching. The main contributions of this paper are summarized below:(i)
We propose an aerial-ground integrated network architecture in disaster emergency areas. The aerial sensing access layer can sense ground information and provide wireless coverage for the ground at the same time(ii)
We design a new mobility model that considers disaster scenarios based on the 3D Gaussian Markov mobility model. Its movement direction integrates the movement information of the disaster-stricken people on the ground in order to collect information more accurately(iii)
We design a routing metric for UAV dynamics and path instability and establish a multipath stability optimization function combined with spectrum sensing(iv)
We design a multidimensional hypergraph matching model for multipath selection and propose a multidimensional hypergraph matching algorithm. Firstly, the aerial sensing access layer is divided into multiple regions by aerial base station (ABS) arranged at equal intervals ask-partite hypergraph. Taking UAVs as hypergraph vertices, the entire path connecting each region and simultaneously channel selection is taken as the hyperedge. The stability routing metric, which taking UAV dynamics, link quality, and distance between nodes into account, is taken as the weight of hyperedge. Finally, the problem of maximum optimization stability is transformed into the problem of maximum weight of multidimensional hypergraph matchingThe rest of this paper is organized as follows: In Section2, we review the related work. In Section 3, system model is introduced for multipath stability routing. In Section 4, we transform the problem and propose a (2M-1) D hypergraph matching algorithm. The results and discussion are presented in Section 5. Finally, the conclusions are given in Section 6.


## 2. Related Work

### 2.1. Multipath Routing in UAV Swarms

Here are some literatures introducing multipath routing in UAV swarms from different perspectives. In [6], a network coding-aware multipath routing protocol was proposed, which was based on the ant colony algorithm and achieved end-to-end reliable transmission by maximizing the coding opportunity. In [7], in order to restore, construct, and select k-disjoint paths, an industry-oriented canonical particle swarm (CPS) optimized data delivery framework was introduced. Multiswarm strategy was used to determine the optimal direction of multipath routing. In [8], a new scheme of 3D transformation routing was proposed, which broke down a large topology change into a series of smooth, continuous microadjustments that minimized packet loss (due to path interruptions). In [9], a Lyapunov-based energy-efficient path diversity method for data transfer was proposed in UAV networks. The method first calculated the delay robustness based on path diversity and derived the trade-off relationship between path diversity and power consumption, and the proposed multipath transmission method provided latency tolerance with acceptable power consumption in data transmission. In [10], the authors proposed a SDN and MQTT hybrid structure for battlefield UAV swarms and also developed a QoS-based multi-path routing framework that computes the multiple disjoint paths from source to destination to improve network performance. In [11], a set of multipath scalable video coding routing protocol was designed. Using scalable video coding (SVC) theory, the AODV-multipath (AODVM) with disconnected nodes was extended to enable routing nodes to obtain information such as link parameters. Although the above studies designed UAV swarm multipath routing at different levels, there was little work on multipath stability. In this paper, a new stable routing metric is designed based on the dynamics of the UAV, link quality, and distance between nodes, and a hypergraph matching framework is used to obtain multipath routing.


### 2.2. Routing Based on Hypergraph

Recently, efforts have been made to design routing using hypergraph. In [12], mining algorithms were used in a P2P environment to improve the efficiency of such methods. The proposed method combined clustering and hypergraph, used ECCLAT to build approximate clusters, and found meaningful clusters with slight overlap; this cluster set improved the robustness and availability of query, routing mechanisms in P2P networks extensibility. In [13], it considered the coordination problem of a robot moving on a network. The problem could be solved by introducing social laws, which could be derived from properly defined routing of graphs under the network. Here, it studied the complexity of routing, providing an efficient representation of 2-routable graphs, establishing a correspondence between hypergraph coloring and graph routing. Although there is some work on clustering and coloring applied to routing, it is not suitable for establishing multipath routing.More recently, hypergraph matching theory has begun to be used to design routing. In [14], the authors investigated the shortest path problem in hypergraphs and proposed two algorithms for finding shortest hyperpaths in dynamic networks with dynamic weights and topology. In [15], it studied the problem of finding the optimal path in a single-source single-target accumulative multihop network. By modeling the routing problem in accumulative multihop network as a hypergraph, it deduced a new set of sufficient conditions for optimality, proving useless of the optimality of traditional multihop networks. Although there were some studies on designing routes using hypergraph matching algorithms, few studies had applied hypergraph matching algorithms to multipath routing. We creatively propose a 2M−1-dimensional hypergraph matching algorithm, which is different from the traditional 3D hypergraph matching algorithm.


## 2.1. Multipath Routing in UAV Swarms

Here are some literatures introducing multipath routing in UAV swarms from different perspectives. In [6], a network coding-aware multipath routing protocol was proposed, which was based on the ant colony algorithm and achieved end-to-end reliable transmission by maximizing the coding opportunity. In [7], in order to restore, construct, and select k-disjoint paths, an industry-oriented canonical particle swarm (CPS) optimized data delivery framework was introduced. Multiswarm strategy was used to determine the optimal direction of multipath routing. In [8], a new scheme of 3D transformation routing was proposed, which broke down a large topology change into a series of smooth, continuous microadjustments that minimized packet loss (due to path interruptions). In [9], a Lyapunov-based energy-efficient path diversity method for data transfer was proposed in UAV networks. The method first calculated the delay robustness based on path diversity and derived the trade-off relationship between path diversity and power consumption, and the proposed multipath transmission method provided latency tolerance with acceptable power consumption in data transmission. In [10], the authors proposed a SDN and MQTT hybrid structure for battlefield UAV swarms and also developed a QoS-based multi-path routing framework that computes the multiple disjoint paths from source to destination to improve network performance. In [11], a set of multipath scalable video coding routing protocol was designed. Using scalable video coding (SVC) theory, the AODV-multipath (AODVM) with disconnected nodes was extended to enable routing nodes to obtain information such as link parameters. Although the above studies designed UAV swarm multipath routing at different levels, there was little work on multipath stability. In this paper, a new stable routing metric is designed based on the dynamics of the UAV, link quality, and distance between nodes, and a hypergraph matching framework is used to obtain multipath routing.


## 2.2. Routing Based on Hypergraph

Recently, efforts have been made to design routing using hypergraph. In [12], mining algorithms were used in a P2P environment to improve the efficiency of such methods. The proposed method combined clustering and hypergraph, used ECCLAT to build approximate clusters, and found meaningful clusters with slight overlap; this cluster set improved the robustness and availability of query, routing mechanisms in P2P networks extensibility. In [13], it considered the coordination problem of a robot moving on a network. The problem could be solved by introducing social laws, which could be derived from properly defined routing of graphs under the network. Here, it studied the complexity of routing, providing an efficient representation of 2-routable graphs, establishing a correspondence between hypergraph coloring and graph routing. Although there is some work on clustering and coloring applied to routing, it is not suitable for establishing multipath routing.More recently, hypergraph matching theory has begun to be used to design routing. In [14], the authors investigated the shortest path problem in hypergraphs and proposed two algorithms for finding shortest hyperpaths in dynamic networks with dynamic weights and topology. In [15], it studied the problem of finding the optimal path in a single-source single-target accumulative multihop network. By modeling the routing problem in accumulative multihop network as a hypergraph, it deduced a new set of sufficient conditions for optimality, proving useless of the optimality of traditional multihop networks. Although there were some studies on designing routes using hypergraph matching algorithms, few studies had applied hypergraph matching algorithms to multipath routing. We creatively propose a 2M−1-dimensional hypergraph matching algorithm, which is different from the traditional 3D hypergraph matching algorithm.


## 3. System Model and Problem Formulation

### 3.1. Network Architecture

We consider an aerial-ground integrated distributed underlay cognitive UAV swarm network architecture for ground emergency communications, as shown in Figure1. The air-ground integrated architecture consists of two parts, which are terrestrial layer and aerial sensing access layer.
(i)
In the terrestrial layer, there are still a few ground base stations (GBSs) on the affected ground, which serve as the terrestrial cellular primary network. The primary users (PUs) communicate with base station through licensed orthogonal subchannels. In the case that the GBS does not work properly, the primary users connect with the temporary base station provided by the UAV(ii)
In the aerial sensing access layer, a number of UAVs equipped with radio transceiver, camera function, and GPS modules. UAV swarms are arranged in the air in the form of ad hoc networks. There are two functional categories of UAVs in the aerial that play different roles. One type serves as ABS to provide users with wireless coverage and communication services, and the other type provides ground search through the camera function, and it transmits the sensing information to the ground base station in the form of routingFigure 1
Illustration of aerial-ground integrated architecture for emergency communications.In addition, we consider that the cognitive radio technology of the underlying access is applied in the aerial UAV swarm, and the radio transceiver module can be modulated to any channel that is not occupied by the primary users. In this setting, we study multipath stability routing under cognitive UAV swarm.


### 3.2. Mobility Model

In order to facilitate the representation, we establish a 3D Cartesian coordinate system. We divide the entire aerial layer areaQIl,Iw into M spherical regions, and the region m is defined as Qm, m∈M=1,2,⋯,M. The region Qm is spherical with radius R=Il/2M, and the center of the region is an ABS. Nodes in the UAV swarm are randomly distributed around each ABS. The set of UAVs in the region Qm is represented as Nm=1,2,⋯,Nm. The divided planform regions are shown in Figure 2.Figure 2
Illustration of spherical region division of sensing access layer.We consider modeling the movement of the UAV within each region. In order to more accurately serve the disaster-affected victims on the ground, here, we consider a model that follows the movement of victims on the ground. In fact, in the disaster scenario, the moving direction and path of the victims are difficult to model. But we can roughly divide the victim’s behavior into the following categories [16]:
(i)
In order to prevent secondary injury, the victims may move to safe places. We consider the coordinate of the safe placem for each region to be ssm=xsm,ysm, m∈M, and safe places are represented by a set as S=1,2,⋯,S(ii)
Victims run to the coverage area of the ABS looking for aid. We roughly use the location of the ABS to indicate the coordinate of the victims. ABSs are distributed at the same heightz0. The location of ABS m is denoted as wm=xm,ym,z0, m∈MAfter introducing the victim’s behaviors, we define the UAV mobility model related to the victim’s behavior. To represent the positions of the UAV nodes, we divide the UAV’s flight durationT into K time slots, and the time slots K are large enough so that we can consider the position to be roughly constant at each time slot. The position of UAV at time slot k is denoted as qnmk=xnk,ynk,znk, n∈Nm. For the establishment of the speed model, we refer to the 3D-Gauss-Markov (GM) mobility model [17]. The flight speed of the UAV n at time slot k is expressed as
(1)vnk=βvnk−1+1−βv¯+1−β2Gvk−1,where β∈0,1 is the tuning parameter for changing the degree of randomness. v¯ is the average speed of the UAV when k⟶∞. Gvk−1 is a random variable with a Gaussian distribution.We set the flying direction of the UAVn as the model related to the movement behavior of the victim. The UAV n flying direction and pitch variables represent the actual angles at which the UAV is moving. The direction and vertical pitch at time slot k are denoted as
(2)φnk=φ0+1−βφ¯+1−β2Gφk−1,ϑnk=βϑnk−1+1−βϑ¯+1−β2Gϑk−1,where Gφk−1 and Gϑk−1 are two random variables that follow the independent and unrelated Gaussian distributions. In particular, the φ0 represents the direction of the victim’s movement; specifically, we define the direction selection strategy, and it can be expressed as
(3)φ0=φs=cos−1xsm−xmssm−wmpsφA=cos−1xmwmpAφ0p0,where φ0 selects different directions with different probabilities, respectively. The vector direction φA of the safe place is selected with probability ps, s∈S. The direction of the ABS φA is selected with probability pA, and the direction φ0 is selected with the probability of p0. For ease of understanding, we project the ABS and UAVs on the ground area, and the moving direction selection strategy is shown in Figure 3.Figure 3
Illustration of the moving direction selection strategy.Further, we can obtain the coordinates at time slotk, and it can be expressed by the formula as
(4)xnk=xnk−1+vnk−1cosφnk−1cosϑnk−1,ynk=ynk−1+vnk−1sinφnk−1cosϑnk−1,znk=znk−1+vnk−1sinϑnk−1.


### 3.3. Stability Model

According to the movement model, the distance between nodes and the speed of the UAV can be easily obtained. In order to obtain stable multipath routing, according to the dynamic characteristics of the UAV, we jointly consider the relative speed and distance between nodes, as well as the link quality.In order to represent the link quality, we introduce the received signal strength indicator (RSSI) to indicate the quality of the link, because compared with other link quality, RSSI can judge the link quality more quickly and accurately [18]. The link quality between the UAV i in the region m and the UAV j in the adjacent next region is specifically expressed as
(5)Ii,jm=∑x=1NRSRSx2NRS−∑x=1NRSRSxNRS2,where NRS represents the total number of RSSI samples, and RSx represents the value of the xth sample of RSSI.The dynamics of UAVi and UAV j within region m is expressed as a velocity-dependent exponential function
(6)Vi,jm=1−ηevi−vj/v¯,where vi is the node speed in the region m and vj is the node speed in the next region, v¯ indicates the average speed of the UAV, and η∈0,1 is a tuning parameter, ensuring that Vi,jm is a positive number.Combining the indicators defined above, we further define the degree of stability metric as follows(7)Ui,jmk=μIi,jm+λVi,jmςdi,jk,whereμ and λ are weighting factors, and μ+λ=1, in order to balance the proportion of dynamic and link quality between UAVs. dijk represents the distance between UAVs i and j at time slot k. Specially, ς is the distance factor to adjust the whole formula, so that the obtained ratio is in the normal value range. According to the mobility model, we can express the distance between UAVs i and j at time slot k as
(8)di,jk=qik−qjkxik−xjk2+yik−yjk2+zik−zjk2.We use a binary variable to represent the relationship between nodei and node j at time slot k in the mth region
(9)ai,jmk=1,ifUAViassociatedwithUAVj0,otherwise.


### 3.4. Spectrum Sensing Model

After designing the selection metric between nodes, channel selection is also required during data transmission. For the primary network on the ground, we divide the entire RB intoCm authorized subchannels in the mth region and set as Cm=1,2,⋯,Cm. According to the ON/OFF model [19], the ON state means the channel is occupied, the OFF state means the channel is idle, and the CU can access only when the authorized channel is OFF state.Here, we set the prior probability that the channelc is idle at time slot k as ρcmk. In practice, the available probability of each channel is expressed as
(10)ρcmk⟵ρcmk+δc,δc=12ρcmk−1−ρcmk−2+ρcmk−2−ρcmk−3,where δc is the correction factor for the channel availability probability.We use binary variables to represent the relationship between nodei and channel c at time slot k in the mth region
(11)bi,cmk=1,&ifUAViassociatedwithchannelc0,otherwise.Suppose we buildL paths, the multipath stability from the source node to the destination node at time slot k is expressed as
(12)Sk=∑m=1M∑i=1L∑j=1Nmai,jmkUijmk+∑c=1Cmbi,cmkρcmk.


### 3.5. Problem Formulation

After the above settings, we aim to maximize paths stability over timeT by optimizing the binary variables of M hop UAV nodes and M−1 channels selection. This optimization problem can be formulated mathematically as
(13)maxai,jmk,bi,cmk∑k=1KSk,(13a)s.t.∑i=1Lai,jmk≤L,∀j,k,(13b)∑j=1Nmai,jmk≤1,∀i,k,(13c)∑m=1Mai,jmk≤M,∀i,j,(13d)∑i=1Lbi,cmk≤L,∀c,k,(13e)∑c=1Cmbi,cmk≤1,∀i,k,(13f)∑m=1Mbi,cmk≤M,∀i,c,(13g)ai,jmk,bi,cmk∈0,1,∀i,j,c,m,where (13a) and (13b) are UAV selection relationship constraints between adjacent regions. (13d) and (13e) are channel selection relationship constraints between adjacent regions. And (13c) and (13f) are the constraints of the number of the regions. Finally, (13g) reflects the constraints on the binary variables of UAV selection and channel selection.


## 3.1. Network Architecture

We consider an aerial-ground integrated distributed underlay cognitive UAV swarm network architecture for ground emergency communications, as shown in Figure1. The air-ground integrated architecture consists of two parts, which are terrestrial layer and aerial sensing access layer.
(i)
In the terrestrial layer, there are still a few ground base stations (GBSs) on the affected ground, which serve as the terrestrial cellular primary network. The primary users (PUs) communicate with base station through licensed orthogonal subchannels. In the case that the GBS does not work properly, the primary users connect with the temporary base station provided by the UAV(ii)
In the aerial sensing access layer, a number of UAVs equipped with radio transceiver, camera function, and GPS modules. UAV swarms are arranged in the air in the form of ad hoc networks. There are two functional categories of UAVs in the aerial that play different roles. One type serves as ABS to provide users with wireless coverage and communication services, and the other type provides ground search through the camera function, and it transmits the sensing information to the ground base station in the form of routingFigure 1
Illustration of aerial-ground integrated architecture for emergency communications.In addition, we consider that the cognitive radio technology of the underlying access is applied in the aerial UAV swarm, and the radio transceiver module can be modulated to any channel that is not occupied by the primary users. In this setting, we study multipath stability routing under cognitive UAV swarm.


## 3.2. Mobility Model

In order to facilitate the representation, we establish a 3D Cartesian coordinate system. We divide the entire aerial layer areaQIl,Iw into M spherical regions, and the region m is defined as Qm, m∈M=1,2,⋯,M. The region Qm is spherical with radius R=Il/2M, and the center of the region is an ABS. Nodes in the UAV swarm are randomly distributed around each ABS. The set of UAVs in the region Qm is represented as Nm=1,2,⋯,Nm. The divided planform regions are shown in Figure 2.Figure 2
Illustration of spherical region division of sensing access layer.We consider modeling the movement of the UAV within each region. In order to more accurately serve the disaster-affected victims on the ground, here, we consider a model that follows the movement of victims on the ground. In fact, in the disaster scenario, the moving direction and path of the victims are difficult to model. But we can roughly divide the victim’s behavior into the following categories [16]:
(i)
In order to prevent secondary injury, the victims may move to safe places. We consider the coordinate of the safe placem for each region to be ssm=xsm,ysm, m∈M, and safe places are represented by a set as S=1,2,⋯,S(ii)
Victims run to the coverage area of the ABS looking for aid. We roughly use the location of the ABS to indicate the coordinate of the victims. ABSs are distributed at the same heightz0. The location of ABS m is denoted as wm=xm,ym,z0, m∈MAfter introducing the victim’s behaviors, we define the UAV mobility model related to the victim’s behavior. To represent the positions of the UAV nodes, we divide the UAV’s flight durationT into K time slots, and the time slots K are large enough so that we can consider the position to be roughly constant at each time slot. The position of UAV at time slot k is denoted as qnmk=xnk,ynk,znk, n∈Nm. For the establishment of the speed model, we refer to the 3D-Gauss-Markov (GM) mobility model [17]. The flight speed of the UAV n at time slot k is expressed as
(1)vnk=βvnk−1+1−βv¯+1−β2Gvk−1,where β∈0,1 is the tuning parameter for changing the degree of randomness. v¯ is the average speed of the UAV when k⟶∞. Gvk−1 is a random variable with a Gaussian distribution.We set the flying direction of the UAVn as the model related to the movement behavior of the victim. The UAV n flying direction and pitch variables represent the actual angles at which the UAV is moving. The direction and vertical pitch at time slot k are denoted as
(2)φnk=φ0+1−βφ¯+1−β2Gφk−1,ϑnk=βϑnk−1+1−βϑ¯+1−β2Gϑk−1,where Gφk−1 and Gϑk−1 are two random variables that follow the independent and unrelated Gaussian distributions. In particular, the φ0 represents the direction of the victim’s movement; specifically, we define the direction selection strategy, and it can be expressed as
(3)φ0=φs=cos−1xsm−xmssm−wmpsφA=cos−1xmwmpAφ0p0,where φ0 selects different directions with different probabilities, respectively. The vector direction φA of the safe place is selected with probability ps, s∈S. The direction of the ABS φA is selected with probability pA, and the direction φ0 is selected with the probability of p0. For ease of understanding, we project the ABS and UAVs on the ground area, and the moving direction selection strategy is shown in Figure 3.Figure 3
Illustration of the moving direction selection strategy.Further, we can obtain the coordinates at time slotk, and it can be expressed by the formula as
(4)xnk=xnk−1+vnk−1cosφnk−1cosϑnk−1,ynk=ynk−1+vnk−1sinφnk−1cosϑnk−1,znk=znk−1+vnk−1sinϑnk−1.


## 3.3. Stability Model

According to the movement model, the distance between nodes and the speed of the UAV can be easily obtained. In order to obtain stable multipath routing, according to the dynamic characteristics of the UAV, we jointly consider the relative speed and distance between nodes, as well as the link quality.In order to represent the link quality, we introduce the received signal strength indicator (RSSI) to indicate the quality of the link, because compared with other link quality, RSSI can judge the link quality more quickly and accurately [18]. The link quality between the UAV i in the region m and the UAV j in the adjacent next region is specifically expressed as
(5)Ii,jm=∑x=1NRSRSx2NRS−∑x=1NRSRSxNRS2,where NRS represents the total number of RSSI samples, and RSx represents the value of the xth sample of RSSI.The dynamics of UAVi and UAV j within region m is expressed as a velocity-dependent exponential function
(6)Vi,jm=1−ηevi−vj/v¯,where vi is the node speed in the region m and vj is the node speed in the next region, v¯ indicates the average speed of the UAV, and η∈0,1 is a tuning parameter, ensuring that Vi,jm is a positive number.Combining the indicators defined above, we further define the degree of stability metric as follows(7)Ui,jmk=μIi,jm+λVi,jmςdi,jk,whereμ and λ are weighting factors, and μ+λ=1, in order to balance the proportion of dynamic and link quality between UAVs. dijk represents the distance between UAVs i and j at time slot k. Specially, ς is the distance factor to adjust the whole formula, so that the obtained ratio is in the normal value range. According to the mobility model, we can express the distance between UAVs i and j at time slot k as
(8)di,jk=qik−qjkxik−xjk2+yik−yjk2+zik−zjk2.We use a binary variable to represent the relationship between nodei and node j at time slot k in the mth region
(9)ai,jmk=1,ifUAViassociatedwithUAVj0,otherwise.


## 3.4. Spectrum Sensing Model

After designing the selection metric between nodes, channel selection is also required during data transmission. For the primary network on the ground, we divide the entire RB intoCm authorized subchannels in the mth region and set as Cm=1,2,⋯,Cm. According to the ON/OFF model [19], the ON state means the channel is occupied, the OFF state means the channel is idle, and the CU can access only when the authorized channel is OFF state.Here, we set the prior probability that the channelc is idle at time slot k as ρcmk. In practice, the available probability of each channel is expressed as
(10)ρcmk⟵ρcmk+δc,δc=12ρcmk−1−ρcmk−2+ρcmk−2−ρcmk−3,where δc is the correction factor for the channel availability probability.We use binary variables to represent the relationship between nodei and channel c at time slot k in the mth region
(11)bi,cmk=1,&ifUAViassociatedwithchannelc0,otherwise.Suppose we buildL paths, the multipath stability from the source node to the destination node at time slot k is expressed as
(12)Sk=∑m=1M∑i=1L∑j=1Nmai,jmkUijmk+∑c=1Cmbi,cmkρcmk.


## 3.5. Problem Formulation

After the above settings, we aim to maximize paths stability over timeT by optimizing the binary variables of M hop UAV nodes and M−1 channels selection. This optimization problem can be formulated mathematically as
(13)maxai,jmk,bi,cmk∑k=1KSk,(13a)s.t.∑i=1Lai,jmk≤L,∀j,k,(13b)∑j=1Nmai,jmk≤1,∀i,k,(13c)∑m=1Mai,jmk≤M,∀i,j,(13d)∑i=1Lbi,cmk≤L,∀c,k,(13e)∑c=1Cmbi,cmk≤1,∀i,k,(13f)∑m=1Mbi,cmk≤M,∀i,c,(13g)ai,jmk,bi,cmk∈0,1,∀i,j,c,m,where (13a) and (13b) are UAV selection relationship constraints between adjacent regions. (13d) and (13e) are channel selection relationship constraints between adjacent regions. And (13c) and (13f) are the constraints of the number of the regions. Finally, (13g) reflects the constraints on the binary variables of UAV selection and channel selection.


## 4. Problem Transformation and Algorithm Design

The problem in (13) is a 0-1 integer linear programming problem, and its goal is to find a binary variable assigned a value of 0 or 1 to maximize the objective function and satisfy all constraints. Such problems are NP-complete. When the number of channels and the number of UAVs in the system increase, it will be very difficult to solve by traditional methods. Therefore, this section transforms the problem of M UAV nodes and M−1 channels selecting in M regions in each path into a weighted 2M−1-uniform hypergraph model.


### 4.1. Hypergraph Construction

Definition 1 (Hypergraph).
A hypergraph can be defined asH=V,E, where V is the set of vertices and E is the set of hyperedges. The weighted hypergraph is expressed as H=V,E,w, where the weight of each hyperedge is denoted as w=Ue, for e=1,2,⋯,E
In the aerial layer, we divideM regions, and the set of nodes in region m is represented as Nm. The set of all UAVs used to sense the ground in the aerial layer is represented as the set of hypergraph vertices; i.e., V=∑m=1MNm. Since the stability of nodes and the available probabilities of channels are dynamically changing in different time slots, we also introduce time slots into the hypergraph. Then, we define the weight of the hyperedge as the stability metric; i.e., Uek=Ui,jmk+ρcmk.Definition 2 (J-Partite Hypergraph).
A hypergraph is called aJ-partite hypergraph if its vertices can be divided into J subsets. Therein, the hyperedge contains a vertex selected from each subset. It is also called J-dimensional hypergraph.In our hypergraph, we have divided the UAVs intoM regions. In addition, we need to make M−1 channels selection, and then, it actually divides into 2M−1 subsets of vertices. Therefore, our optimization problem can be viewed as a 2M−1-dimensional hypergraph matching problem.


### 4.2. Algorithm Design

For2M−1-dimensional hypergraph matching, finding the largest subset of weights for vertex disjoint hyperedges is a NP-hard problem. Inspired by the idea of local search with increasing approximation ratio in polynomial time [20], and in this section, we adopt local search to design hypergraph matching algorithms to find suboptimal solutions.Definition 3 (Conflict Graph).
We denote the conflict graph asD=G,Z, which is a special graph transformed from hypergraph. Anyone’s vertices g∈G correspond to one hyperedge e∈E. If there is at least one common vertex in the hyperedges e1 and e2, and then, g1 and g2 correspond to the same situation, which can be expressed as Ng1=g2.As shown in Figure4, we show the selection of channels and nodes between region m and region m+1. To illustrate the problem, we only take four nodes and four channels in per region. Here, we define the original hypergraph as the graph connecting all vertices to form hyperedges. According to the constraints, a feasible hypergraph can be obtained. Further, initialize a set of independent hyperedge sets as GA in the conflict graph from the feasible hypergraph. We define GB=G−GA in the conflict graph as the set that contain all adjacent vertices in the set of GA [21].Definition 4 (ω−claw).
Theω-claw graph is an induced subgraph of the conflict graph, denoted as Lω, consisting of ω independent vertices, called as talons [22], which forms an independent set of YLω. The center vertex is in set GA and the ω independent claw is in set GB.
In the2M−1-dimensional matching problem, the central vertex of the conflict graph is adjacent to at most 2M−1 independent vertices. There are 2M−1 subsets in total, and each subset can be repeated once so that the adjacent vertices can be obtained, which means ω≤2M−1.Figure 4
Illustration of original hypergraph, feasible hypergraph and conflict graph.To obtain the initial and independent setGA, the Algorithm 1 uses greedy algorithm to get the set of disjoint hyperedges with the largest weight. Then, we iteratively search for ω-claw by local search, a 2M−1-dimensional hypergraph matching algorithm is proposed. In Algorithm 2, the GA obtained in Algorithm 1 is firstly sorted in ascending order, and the search starts from the vertex with the smallest weight. Since Definition 4 analyzes ω≤2M−1, the algorithm starts to search from the case of ω=1, in the set GB finding vertices adjacent to the first vertex. We use NYLω,GAk to represent the adjacent independent set YLω with ω vertices in an initial independent set GAk [20]. If ϒGAk−NYLω,GAk∪YLω≻ϒGAk means that there is better stability in the adjacent vertices, then update GA. We can search up to the case where ω=2M−1. According to the output of Algorithm 2, we can obtain the nodes disjoint largest subset of weights. Besides, the maximum stability at time slot k is expressed as YGAk. Therefore, the stability of the system during the entire time T can be expressed as ∑k=1KYGAk. Compared with the traditional greedy algorithm, the proposed algorithm in this paper adopts the local search method and improves the overall search performance of the algorithm by identifying the ω-claw graph.Algorithm 1: Greedy Algorithm.
Input: Stability UekOutput: The selected initialization set GAk1: Transform the original hypergraph into a feasible hypergraph.2: Transform the feasible hypergraph into a conflict graph.3: Generate conflict graphD=G,Z with vertex set Z based on H=V,E.4: SetGAk=∅.5:whileGA≠Ldo6: Select a vertexz∗∈Z with maximum-weight of vertex in D=G,Z.7: SetNz∗k=0. Find all the adjacent vertices of z∗. And let Nz∗k be the set of adjacent vertices of z∗.8: SetGAk=GAk∪z∗ and Z=Z−z∗∪Nz∗k9:end whileAlgorithm 2:2M−1-Dimensional Hypergraph Matching Algorithm.
Input: Stability UekOutput: Hypergraph matching E∗1: Transform the original hypergraph into a feasible hypergraph.2: Transform the feasible hypergraph into a conflict graph.3: Use Algorithm5 to get the initial independent set GAk.4: Sort all the vertices ofGAk in an ascending order based on the weight of every vertex. Set i=1.5:fork=1 to Kdo6:whileω≤2M−1do7:  Search forω-claw Lω from the set D.8:if there exists ω-claw YLω satisfying9:ϒGAk−NYLω,GAk∪YLω≻ϒGAkthen10:GAk⟵GAk−NYLω,GAk∪YLω, go to step 4.11:else12:ω=ω+113:end if14:end while15:ifi≤GAthen16:i=i+1.17:end if18:end for


### 4.3. Complexity Analysis

The computational complexity of Algorithm2 mainly comes from the complexity of the inner layer Algorithm 1 and the cyclic computational complexity of the outer layer hypergraph matching algorithm. First, we analyze the computational complexity of Algorithm 1. We need to find L independent vertices. From Algorithm 1, we can select a hyperedge with the largest weight in Z each time. In the worst case, the computational complexity is OLZ. In Algorithm 2, the computational complexity of line 3 is the computational complexity of Algorithm 1. The line 4 is sorted in ascending order, and the computational complexity is OMlogM. The line 5 is a for loop, and the inner layer seeks ω-claw, and we set the number of adjacent vertices to GB, and then, the complexity is OKGB2M−1 [23]. Therefore, the computational complexity of the entire algorithm is OLZ+MlogM+KGB2M−1, which can be approximated as OKGB2M−1.


## 4.1. Hypergraph Construction

Definition 1 (Hypergraph).
A hypergraph can be defined asH=V,E, where V is the set of vertices and E is the set of hyperedges. The weighted hypergraph is expressed as H=V,E,w, where the weight of each hyperedge is denoted as w=Ue, for e=1,2,⋯,E
In the aerial layer, we divideM regions, and the set of nodes in region m is represented as Nm. The set of all UAVs used to sense the ground in the aerial layer is represented as the set of hypergraph vertices; i.e., V=∑m=1MNm. Since the stability of nodes and the available probabilities of channels are dynamically changing in different time slots, we also introduce time slots into the hypergraph. Then, we define the weight of the hyperedge as the stability metric; i.e., Uek=Ui,jmk+ρcmk.Definition 2 (J-Partite Hypergraph).
A hypergraph is called aJ-partite hypergraph if its vertices can be divided into J subsets. Therein, the hyperedge contains a vertex selected from each subset. It is also called J-dimensional hypergraph.In our hypergraph, we have divided the UAVs intoM regions. In addition, we need to make M−1 channels selection, and then, it actually divides into 2M−1 subsets of vertices. Therefore, our optimization problem can be viewed as a 2M−1-dimensional hypergraph matching problem.


## 4.2. Algorithm Design

For2M−1-dimensional hypergraph matching, finding the largest subset of weights for vertex disjoint hyperedges is a NP-hard problem. Inspired by the idea of local search with increasing approximation ratio in polynomial time [20], and in this section, we adopt local search to design hypergraph matching algorithms to find suboptimal solutions.Definition 3 (Conflict Graph).
We denote the conflict graph asD=G,Z, which is a special graph transformed from hypergraph. Anyone’s vertices g∈G correspond to one hyperedge e∈E. If there is at least one common vertex in the hyperedges e1 and e2, and then, g1 and g2 correspond to the same situation, which can be expressed as Ng1=g2.As shown in Figure4, we show the selection of channels and nodes between region m and region m+1. To illustrate the problem, we only take four nodes and four channels in per region. Here, we define the original hypergraph as the graph connecting all vertices to form hyperedges. According to the constraints, a feasible hypergraph can be obtained. Further, initialize a set of independent hyperedge sets as GA in the conflict graph from the feasible hypergraph. We define GB=G−GA in the conflict graph as the set that contain all adjacent vertices in the set of GA [21].Definition 4 (ω−claw).
Theω-claw graph is an induced subgraph of the conflict graph, denoted as Lω, consisting of ω independent vertices, called as talons [22], which forms an independent set of YLω. The center vertex is in set GA and the ω independent claw is in set GB.
In the2M−1-dimensional matching problem, the central vertex of the conflict graph is adjacent to at most 2M−1 independent vertices. There are 2M−1 subsets in total, and each subset can be repeated once so that the adjacent vertices can be obtained, which means ω≤2M−1.Figure 4
Illustration of original hypergraph, feasible hypergraph and conflict graph.To obtain the initial and independent setGA, the Algorithm 1 uses greedy algorithm to get the set of disjoint hyperedges with the largest weight. Then, we iteratively search for ω-claw by local search, a 2M−1-dimensional hypergraph matching algorithm is proposed. In Algorithm 2, the GA obtained in Algorithm 1 is firstly sorted in ascending order, and the search starts from the vertex with the smallest weight. Since Definition 4 analyzes ω≤2M−1, the algorithm starts to search from the case of ω=1, in the set GB finding vertices adjacent to the first vertex. We use NYLω,GAk to represent the adjacent independent set YLω with ω vertices in an initial independent set GAk [20]. If ϒGAk−NYLω,GAk∪YLω≻ϒGAk means that there is better stability in the adjacent vertices, then update GA. We can search up to the case where ω=2M−1. According to the output of Algorithm 2, we can obtain the nodes disjoint largest subset of weights. Besides, the maximum stability at time slot k is expressed as YGAk. Therefore, the stability of the system during the entire time T can be expressed as ∑k=1KYGAk. Compared with the traditional greedy algorithm, the proposed algorithm in this paper adopts the local search method and improves the overall search performance of the algorithm by identifying the ω-claw graph.Algorithm 1: Greedy Algorithm.
Input: Stability UekOutput: The selected initialization set GAk1: Transform the original hypergraph into a feasible hypergraph.2: Transform the feasible hypergraph into a conflict graph.3: Generate conflict graphD=G,Z with vertex set Z based on H=V,E.4: SetGAk=∅.5:whileGA≠Ldo6: Select a vertexz∗∈Z with maximum-weight of vertex in D=G,Z.7: SetNz∗k=0. Find all the adjacent vertices of z∗. And let Nz∗k be the set of adjacent vertices of z∗.8: SetGAk=GAk∪z∗ and Z=Z−z∗∪Nz∗k9:end whileAlgorithm 2:2M−1-Dimensional Hypergraph Matching Algorithm.
Input: Stability UekOutput: Hypergraph matching E∗1: Transform the original hypergraph into a feasible hypergraph.2: Transform the feasible hypergraph into a conflict graph.3: Use Algorithm5 to get the initial independent set GAk.4: Sort all the vertices ofGAk in an ascending order based on the weight of every vertex. Set i=1.5:fork=1 to Kdo6:whileω≤2M−1do7:  Search forω-claw Lω from the set D.8:if there exists ω-claw YLω satisfying9:ϒGAk−NYLω,GAk∪YLω≻ϒGAkthen10:GAk⟵GAk−NYLω,GAk∪YLω, go to step 4.11:else12:ω=ω+113:end if14:end while15:ifi≤GAthen16:i=i+1.17:end if18:end for


## 4.3. Complexity Analysis

The computational complexity of Algorithm2 mainly comes from the complexity of the inner layer Algorithm 1 and the cyclic computational complexity of the outer layer hypergraph matching algorithm. First, we analyze the computational complexity of Algorithm 1. We need to find L independent vertices. From Algorithm 1, we can select a hyperedge with the largest weight in Z each time. In the worst case, the computational complexity is OLZ. In Algorithm 2, the computational complexity of line 3 is the computational complexity of Algorithm 1. The line 4 is sorted in ascending order, and the computational complexity is OMlogM. The line 5 is a for loop, and the inner layer seeks ω-claw, and we set the number of adjacent vertices to GB, and then, the complexity is OKGB2M−1 [23]. Therefore, the computational complexity of the entire algorithm is OLZ+MlogM+KGB2M−1, which can be approximated as OKGB2M−1.


## 5. Results and discussions

Simulation results are presented in this section to evaluate the proposed multipath stability disjoint routing architecture in disasters. The horizontal projection of the aerial sensing access layer is a two-dimensional rectangular areaQ10km,10km. We divide the whole area into M∈10,20 regions. The number of the UAVs in each region are Nm=30,40, and the flight during time T=12s. The entire RB is divided into C=Cm=50 channels.Regarding the mobility model, each of our region generatesS=2 safe places within the range of the ground. The altitude of ABS is set to z0=50m. The location of the ABS 1 w1=50m,70m,50m. The remaining ABSs are equally spaced in other locations. For the UAV, the tuning parameter β=0.5. The average speed of UAV is set to v¯=12m/s. The direction and the vertical pitch are set to φ¯=1/5πrad and ϑ¯=0. Gvk−1, Gφk−1 and Gϑk−1 fit Gaussian (normal) distribution 0,1, respectively. As for the direction selection strategy, we set the p1=0.2, p2=0.2, pA=0.2, and p0=0.2. φ0 is initialized to zero.As for the stability model, the number of RSSI samples isNRS=50, and the value of each sample lies in the interval −70dBm,−50dBm [18]. The tuning parameter η=0.4. Besides, we set the weighting factors μ=λ=0.5 and the distance factor ς=0.01. Finally, we set the channel prior probability ρcmk∈0.3,0.7.The accumulated stability over the entire time period is denoted as∑k=1KYGAk. In order to express the relationship between the number of time slots and the stability, the average time slot stability is defined as
(14)τ=∑k=1KYGAkK.In order to express the relationship between the number of paths and the stability, the average path stability is defined, which is used to represent the average stability of each path, i.e.,(15)∆=∑k=1KYGAkL.The average regional stability is further defined to express the relationship between the number of regions and the stability, and the average stability within each region is expressed as(16)Ω=∑k=1KYGAkM.In order to better evaluate the performance of the proposed algorithm, this section considers the following two contrasting schemes for comparison: the first is the greedy algorithm (GA), and the second is the proposed hypergraph matching algorithm (HMA) when the number of nodes within each region isNm=40 and Nm=30, respectively.The relationship between the number of slots and the average slot stability is shown in Figure5. When the number of channels C=40, the number of paths and the number of regions are selected L=3 and M=10, respectively. We compared that the average time slot stability of the GA with the number of nodes in each region is Nm=30 and the proposed algorithm HMA with the number of nodes in each region is Nm=40 and Nm=30, respectively. It can be seen from the figure that with the increase of the number of time slots, the average time slot stability of the three curves first shows an increasing trend and gradually tends to remain unchanged. The intuitive explanation of this trend can be as follows: the total time remains unchanged, and as the number of time slots increases, the time slot length becomes shorter, more reliable information can be obtained in time, and the performance of obtaining the overall stability of the path is more accurate. In addition, it can be found that the performance of the proposed HMA is better than that of the GA, because by combining ω-claw, the proposed algorithm will try its best to find a suboptimal solution, which can improve the overall search performance of the system. For the HMA to obtain higher stability at Nm=40, it can be explained that when the number of nodes in the same area increases, the system has more choices to choose nodes with better stability, and then a higher stable path can be obtained. The results show that the system can obtain better performance by appropriately increasing the number of time slots, which highlights the importance of properly adjusting the number of time slots.Figure 5
The average time slot stability versus the number of times slots.In Figure6, the relationship between the number of channels and the accumulated stability is shown. When K=50, M=10, and L=3, we compared that the accumulated stability of the GA with the number of nodes in each region is Nm=30 and the proposed algorithm HMA with the number of nodes in each region is Nm=40 and Nm=30, respectively. It can be seen from the figure that with the increase of the number of channels, the curves of the proposed algorithm HMA and GA show a positive correlation trend, but with the increase of the number of channels, the change of the accumulated stability is not obvious. This is because as the number of channels increases, the number of channels available for nodes to select increases, and nodes can select more channels with high idle probability. However, because the channel availability probability has a low weight in the optimization function, so the effect will not be obvious. The results show that the influence of the number of channels on the system should be properly considered when setting up an aerial UAV swarm.Figure 6
The accumulated stability versus the number of channels.In Figure7, the relationship between the number of paths and the average path stability is shown. When K=50, M=10, and C=40, we compared that the average path stability of the GA with the number of nodes in each region is Nm=30, and the proposed algorithm HMA with the number of nodes in each region is Nm=40 and Nm=30, respectively. It can be seen from the figure that with the increase of the number of paths, the average path stability shows a rapidly decreasing trend when the number of paths is 4. This is because at the beginning, with the increase of the number of paths, the paths will tend to be diversified, and the stability will not be reduced. When the number of paths set exceeds the range, it exceeds the capacity of the system, resulting in the stability of the extra paths, so it starts to decrease, thereby reducing the average path stability. This result shows that although multipath can improve the diversity of paths and improve the stability, it is not the more the better, and too much will have a negative effect. The number of multipaths should be set reasonably according to the capacity of the system, so that the system reaches the optimum.Figure 7
The average path stability versus the number of paths.In Figure8, we show the relationship between the number of regions and the average region stability. When K=50, C=40, and L=3, we compared that the average region stability of the GA with the number of nodes in each region is Nm=30 and the proposed algorithm HMA with the number of nodes in each region is Nm=40 and Nm=30, respectively. It can be found from the figure that with the increase of the number of regions, the average regional stability shows an increasing trend and gradually becomes stable. This is because with the increase of the number of regions, the distance between the regions will become smaller under the premise that the area of the overall disaster remains unchanged. The distance between nodes in the stability index will become smaller as a whole. The distance and stability are inversely proportional. The closer it is, the higher the stability. The results show that the influence of the number of divided regions on the system performance should be reasonably considered when setting the system.Figure 8
The average region stability versus the number of regions.


## 6. Conclusions

This paper studies the multipath routing problem of aerial cognitive UAV swarm in the emergency communication scenario, and we transform the entire multipath stability problem into a 0-1 integer linear programming problem. Considering the channel situation, through the redesign of the entire network structure, the2M−1-dimensional disjoint hypergraph matching problem is established. The suboptimal solution is found by designing a hypergraph matching algorithm using local search. The results show that the algorithm has better performance, and compared with the greedy algorithm, it can obtain multipath with better stability. In the practical communication systems, we have two remarks drawn from both the theoretical analysis and numerical results. A larger number of regions and a reasonably number of paths are beneficial to design the practical cognitive UAV swarm multipath stability routing system.



---

*Source: 6783041-2022-08-01.xml*