Efficient Radio Channel Allocation in Integrated mmWave/ Sub-6GHz UAV-Assisted Disaster Relief Networks

)e unmanned aerial vehicle(UAV-) assisted sub-6GHz disaster relief networks cannot meet high-speed transmission requirements. In this paper, the millimeter wave (mmWave) frequency band is combined with the sub-6GHz frequency band to build a high-speed UAV-assisted disaster relief network. However, the high propagation path loss of mmWave signals usually needs to be compensated by beamforming, where the ground-facing beam of each UAV is the desired receiving beam of ground user information.)e different channels need to be allocated to a single UAV so that this kind of beam can be used simultaneously by different ground users to communicate with this UAV. Also, the other UAVs should reuse these channels as much as possible to save spectrum resources. In this paper, the beamforming training (BFT) mechanism is firstly used to obtain the signal-to-noise ratio (SNR) values of all possible links between ground terminals and UAVs, which are used to estimate these links’ energy efficiency.)en, an interference graph construction algorithm is proposed to identify the links that cannot be used simultaneously in the same channel according to the system energy efficiency. Finally, an iterative channel allocation algorithm is designed to allocate new channels to eliminate the edges of the interference graph, so that the links obtained by the BFTprocess can be used simultaneously as much as possible under the constraint of the number of channels. )e simulation results show that our proposed scheme can achieve the shorter average convergence time, the higher data rate (or the lower data loss rate), and the higher energy efficiency.


Introduction
Over the last decade, various natural disasters have taken place in many countries around the world, causing huge economic losses and even threatening human lives by affecting the emergency information diffusion. Since the base stations (BSs) play an important role in managing resource allocation of wireless communication networks, the communication performance will be degraded seriously if the BSs are damaged due to an unexpected natural disaster [1]. In the absence of communication networks, some crucial tasks (e.g., rapid emergency information diffusion, effective disaster management, and rapid response) cannot be completed in time, and victims and rescue workers (or teams) are unable to communicate with each other or the outside world, which will seriously hinder the timely and effective disaster assistance within 72 hours after the disaster occurred and thus result in greater loss of personnel and property [2]. erefore, in post-disaster areas, it is necessary to repair the damaged BSs on the spot (if possible) or find out alternative methods to recover the emergency communication networks as quickly as possible.
rough serving as flying BSs to offer connection service to ground terminals that are losing connection from ground BSs [3,4], unmanned aerial vehicles (UAVs) can assist the emergency information diffusion since they have strong line-of-sight (LoS) links with ground terminals and the good features such as controllable mobility and low cost. e authors of [5] even focused on UAV-assisted disaster management systems and proposed the network architectures for different types of disasters based on the interaction between UAVs and ground network devices. However, in order to quickly build an effective disaster relief network, reliable channel resources are indispensable. erefore, to make full use of channel resources and reduce channel conflicts [6], it is necessary to manage and allocate them reasonably and efficiently. Also, in order to be fast and robust enough in response to dynamic topology changes, some existing schemes resort to graph theory to solve the dynamic channel allocation problem with tractable complexity.
Generally, a channel allocation problem is transformed into a classical graph coloring problem by assigning appropriate colors (channels) to vertices (users), where any two vertices connected by an edge cannot reuse the same color. In [7], a static and centralized graph coloring approach was introduced into device-to-device (D2D) communications, but it does not apply to highly dynamic graph abstracted from emergency scenarios. By using stochastic learning method, the authors of [8] proposed a fast communicationfree learning (FCFL) algorithm. Although the FCFL can solve the constraint satisfaction problem in a decentralized way, it does not converge fast enough to respond to topology changes. e authors of [9] believed that the main challenges of channel allocation in stochastic and dynamic environments are how to obtain the desired decision results by utilizing the local information and make the convergence speed of channel allocation algorithm keep up with topology changes. erefore, they studied the integration of graph coloring and learning to address the above challenges, where a UAVassisted disaster relief network is modeled as a dynamic topology by using graph theory, and then a dynamic graph coloring algorithm was proposed to address the channel allocation problem of this dynamic topology in a decentralized way.
Although the work in [9] is the first study on UAVassisted emergency communications by introducing the dynamic graph coloring, it only focuses on sub-6 GHz UAVassisted disaster relief networks. e single hop communication in the traditional sub-6 GHz frequency band has long communication distance and good channel quality, but its spectrum resources are very limited. So it is difficult to meet the needs of building a high-rate disaster relief network. Since millimeter wave (mmWave) frequency band has rich spectrum resources to support high data rate wireless transmissions [10][11][12][13], it is a good alternative to sub-6 GHz frequency band.
Due to the high path loss of mmWave communications, the signal attenuation is large, especially when it is blocked [14]. So directional beamforming transmission is helpful to improve the communication quality and enlarge communication range [15]. Beamforming training mechanism [16] is a commonly used and effective method to align beams at both ends of transceivers, which can ensure good communication quality. However, compared with the effective communication distance of the traditional sub-6 GHz frequency band, the mmWave communication distance is still much shorter. According to the estimate in [17], an inter-site distance (ISD) of 75-100 m is required to fully cover independent mmWave deployments. erefore, it is a feasible way to tightly integrate mmWave frequency band with sub-6 GHz frequency band to build a high-rate UAV-assisted disaster relief network. e mmWave band is responsible for ensuring high-speed data transmission while sub-6 GHz band takes charge of providing relatively reliable communication in terms of network control information [18]. e directional beamforming transmission feature in a UAV-assisted mmWave emergency communication network makes each UAV's vertical beam orientation facing the ground become the desired choice of the ground terminals. erefore, the different ground terminals need the different channels to share this desired beam orientation to communicate with this UAV. Also, these channels should be reused by other UAVs as much as possible to improve spectrum efficiency. However, the superposition of mobility and directional beamforming transmission will complicate the channel allocation problem, which needs a new method to cope with interference graph construction and mmWave channel allocation. Our main contributions are given as follows: (1) In order to enable as many terminals on the ground as possible to communicate well with the UAVs simultaneously while minimizing the number of allocated channels, we first use the BFT mechanism in [16] to obtain the SNR values of the communication links between the UAVs and as many ground terminals as possible. en, according to the condition that the system energy efficiency is lower than the given threshold, we propose an interference graph construction algorithm to find out the links that cannot be used concurrently in the same channel. Finally, according to the resulting interference graph, we design an iterative channel allocation algorithm to allocate the reasonable number of channels to the UAVs to reduce the mutual interference, which can ensure that the system energy efficiency is not less than the given threshold.
(2) In order to facilitate mmWave BFT, a training initiation node is needed to coordinate to training process. If the network size is large, multiple initiation nodes can form a hierarchy to deal with it. erefore, the UAVs can act as the training initiation nodes and then the channel allocation problem can be handled with the help of the UAVs. In fact, it is difficult for a pure distributed mode to ensure fast convergence for channel allocation algorithm when the network size gets large. So we adopt a locally centralized mode, which is well integrated with the mmWave BFT process. e rest of this paper is organized as follows. In Section 2, we review the related works in terms of channel allocation problem. e system model and the problem statements are described in Sections 3 and 4, respectively. e design details for the algorithms with respect to interference graph construction and channel allocation are presented in Section 5. We give the convergence analysis in Section 6 and evaluate the simulation results in Section 7. Finally, we conclude this paper in Section 8.

Related Work
To date, numerous channel allocation approaches have been proposed to improve the spectrum efficiency as well as the high-achievable data rate. e authors of [7] utilized hypergraph theory to address the channel allocation problem, which attempts to coordinate the inference between D2D pairs and cellular users. For D2D communication underlying cellular networks, the authors of [19] studied the subchannel allocation scheme based on the simple greedy strategy, while the authors of [20] developed the resource allocation approach to maximize the utility of users. However, those schemes are not suitable for larger-scale networks since they work in a centralized manner. us, some typical decentralized resource allocation schemes [21][22][23] have been developed. In decentralized scenarios, there are two kinds of theories used widely in solving the channel allocation problem: game theory and graph theory. Among the works adopting game theory, the authors of [21] proposed a channel allocation method by using an anticoordination game for interference management in a highly dynamic network. Based on the one-to-one matching theory, the authors of [22] introduced a priority searching algorithm based resource allocation to achieve the utility maximization for D2D communication underlying cellular networks. Combining the matching game and Stackelberg game together, the authors of [23] modeled a decentralized channel-power scheme based on the channel signal information.
Although game theory is a useful tool for solving the resource allocation problem, the computation complexity of the algorithms in the aforementioned works is still high. In particular, the convergence speed and robustness of most game theory methods are not sufficient to deal with the changes of dynamic topology. Hence, in order to cope with this problem, some graph coloring theory-based channel allocation approaches [7,8,24,25] have been proposed. As mentioned above, the work in [7] does not apply to highly dynamic graph, while that in [8] does not converge fast enough to respond to topology changes. To tackle the dynamic spectrum sharing problem, the authors of [24] proposed a dynamic hypergraph coloring approach to suppress the mutual channel interferences at a threshold value. Based on the proposed cooperative D2D communications framework, the authors of [25] adopted graph-based dynamic channel allocation with an attempt to improve the spectral efficiency.
However, the graph coloring theory-based approaches proposed in [7,8,24,25] can only well perform in the case where the BSs can work properly. In the post-disaster areas where BSs are severely damaged, it is necessary to utilize the emerging technologies (e.g., UAV-assisted emergency networks) to restore the communicating services. erefore, some channel allocation methods for UAV-assisted cellular networks emerged [9,[26][27][28]. e authors of [26] aimed to achieve the energy efficiency by jointly optimizing subchannel selection and power control as well as UAV relays deployment, where UAVs just acted as a relay and the mobility of UAVs was not considered. Considering that the movement of UAVs and channel reallocation generate cochannel inferences among different communication links, the authors of [27] proposed channel reallocation methods that dynamically reassign channels to reduce the ripple effect. In [28], the authors introduced nonorthogonal multiple-access to a UAV wireless network architecture and proposed a two-sided matching and swapping algorithm to deal with the problem of limited spectrum resources. Based on local information in decision-making, the authors of [9] studied the integration of graph coloring and learning to adapt to channel allocation under dynamic topology. Although the work in [9] improved the convergence speed, along with those in [7,8,[24][25][26][27][28], it did not apply to channel allocation in integrated mmWave/sub-6 GHz UAV-assisted disaster relief networks.

Network Architecture.
In a typical post-disaster area, most of the ground cellular networks are assumed to go out of service due to the damaged terrestrial BSs. In order to restore the communicating services, we use a set of UAVs to act as flying BSs to send notification messages to survival users and receive the emergency information from them in a point-to-point manner. e cooperation of multiple UAVs can effectively extend the range of network access services. Figure 1 shows the integrated mmWave/sub-6 GHz UAV-assisted disaster relief network architecture adopted in this paper. ere are I slave UAVs (SUAVs), one master UAV (MUAV), and M survival ground stations (GSTAs). We suppose all the SUAVs are within the MUAV's coverage area and controlled by the MUAV through directional mmWave links (i.e., the black beams in Figure 1). Each GSTA chooses one of the UAVs surrounding itself to execute channel access and data transfer through directional mmWave link (i.e., the gray beams in Figure 1). As described in the introduction, directional mmWave links and traditional sub-6 GHz communication links have their own advantages and disadvantages, so we also assume that any GSTA can communicate with any UAV via sub-6 GHz frequency band and all the UAVs communicate with each other via sub-6 GHz frequency band, which can make full use of the advantages of traditional sub-6 GHz frequency band to make up for the shortcomings of mmWave frequency band.
For air-to-ground (A2G), ground-to-air (G2A), and airto-air (A2A) communications in mmWave UAV-assisted disaster relief networks, we consider the commonly used switch-based analog beam pattern [29], which is denoted as follows:

Mobile Information Systems
where φ denotes the beam width of the main lobe in radian, ω represents the beam offset angle to the main lobe in radian, and ε denotes the gain of the side lobe and 0 < ε ≪ 1. For the sake of analysis, we suppose that each of MUAV, SUAV, and GSAT has the limited number of beams at any cross-section, and each beam covers a unique direction in a nonoverlapping manner. us, the number of beams at any crosssection of each of MUAV, SUAV, and GSAT is denoted by n b , which is given by where φ max represents the maximum beam width of MUAV, SUAV, and GSAT. As shown in Figure 1, when a GSAT is covered by multiple UAVs (e.g., GS 3 is covered by SU 1 and SU 2 ), in order to allow the GSAT to choose a suitable UAV (e.g., SU 1 ) among all the UAVs for communication in terms of smaller interference, the GSAT should carry out BFT process with the UAVs around it in advance, where the BFT mechanism in IEEE 802.11ad/ay [30] is a viable option. In addition, based on this BFT mechanism, the authors of [16] proposed an efficient BFT mechanism for dense mmWave networks, which can efficiently establish the directional mmWave links not only between each GSAT and its suitable UAV but also between each SUAV and the MUAV.
Usually, a UAV acting as a flying BS is equipped with multiple radio frequency (RF) chains [31], and its beam facing the ground vertically is the most desired beam of the GSTAs within its coverage. If these GSTAs are assigned to different channels, they transmit data to the same UAV concurrently in this desired beam without interfering with each other.

Time-Slotted Scheme.
In this paper, we discretize the time period Λ into T intervals and represent the set of cumulative intervals as Λ � 1, . . . , t, . . . , T { }. In a three-dimensional Cartesian coordinate model, the real-time coordinate of any GSTA (e.g., GS m ) is represented as During each interval t, we assume that each UAV's position does not change or the change can be ignored. is assumption is feasible and makes sense when the UAVs' flight path can be planned in advance. Even in the environment where the network topology changes randomly, this assumption can be satisfied by adjusting the interval size.
At the beginning of each interval t, the BFT process is started by the MUAV, and then each UAV can get a candidate set of GSTAs based on the BFT results. In this candidate set of GSTAs, the beam performance of each directional mmWave link between the UAV and each GSTA is regarded as energy efficiency of this link, where interference between links is not considered according to the characteristics of BFT mechanism in [16].
If these directional mmWave links are used simultaneously, some links have little or no interference with each other and thus they can reuse the same channel. However, for the links that cannot reuse the same channel, different channels have to be allocated for simultaneous use of these links. e available channel set is denoted by C � c { } C c�1 and the unit bandwidth of one mmWave channel is denoted by B.
Due to the dynamic nature of scenarios, the trained mmWave links may be interrupted. us, the BFT process needs to be redone to find other suitable links. Within each interval, the BFT process, the construction of interference graph, and the mmWave channel allocation should be completed as quickly as possible, so that more part of each interval can be used for data transmission. If the network topology is more dynamic, the interval size should be set smaller, but the channel allocation scheme needs to be faster to keep up with the topology changes, which makes the scheme design more challenging.

Radio Signal Propagation Model and System Energy Efficiency Estimation.
We adopt the following mmWave signal propagation model: In (3), p s i,m is the transmission power at the mmWave link i ⟶ m, while p r i,m is the received power at the mmWave link i ⟶ m; g s i,m is the directional transmitting gain at the mmWave link i ⟶ m, while g r i,m is the directional receiving gain at the mmWave link i ⟶ m; g c i,m is the channel gain at the mmWave link i ⟶ m. When the beam between SU i and GS m is aligned, g s i,m and g r i,m will be calculated by the following formulas [29]: (4b) In (4a) and (4b) φ s i,m represents the transmitter's beam width, while φ r i,m represents the receiver's beam width, and ε represents the gain of the side lobe and 0 < ε ≪ 1. In addition, φ s i,m and φ r i,m represent the main lobe in radian. e channel gain g c i,m is calculated by the following formula [32]: In (5), δ(·) represents the Dirac delta function; τ i,m and χ c i,m represent the propagation delay and the amplitude of the mmWave link i ⟶ m, respectively. τ i,m is calculated by the following formula: In (6), d i,m is the distance of the mmWave link i ⟶ m and c represents the speed of light. When the mmWave link i ⟶ m is a line-of-sight (LOS) link, the amplitude is calculated by the following formula [32]: In (7), λ represents the wavelength, λ � (c/f c ) and f c represents the carrier frequency. When the mmWave link i ⟶ m is a non-line-of-sight (NLOS) link, the amplitude is related to both reflection coefficients and path losses. Due to very high reflection loss in mmWave frequency band [10], only one reflection of a given path is considered, which is calculated by the following formula [32]: In (8), Γ represents the mmWave reflection path's reflection coefficient. e real-time distance of the mmWave link i ⟶ m at each interval t is denoted as d i,m (t), which is estimated by the following formula: e LOS probability of the mmWave link i ⟶ m at each interval t is denoted as P LOS i,m (t), which is estimated by the following formula [9]: In (10), b 1 and b 2 are environment constants, and θ i,m (t) is the elevation angle of the mmWave link i ⟶ m at each interval t, which is estimated by the following formula: e NLOS probability of the mmWave link i ⟶ m at each interval t is denoted by P nlos i,m (t), which is estimated by the following formula [9]: e average channel gain of the mmWave link i ⟶ m at each interval t is denoted as g c i,m (t), which is estimated by the following formula: Based on (13), the received power of the mmWave link i ⟶ m at each interval t is denoted by p r i,m (t), which is estimated by the following formula: For each the mmWave link (e.g., i ⟶ m) determined by the BFT mechanism at each interval t, its SNR is also determined, which is estimated by the following formula: In (15), B and N 0 are the bandwidth of mmWave channel and the power spectral density of background noise, respectively. However, it is not always an effective method to directly measure mmWave channel quality by adopting SNR value. us, considering the pursuit of energy efficiency for high-rate wireless networks, the data rate per unit of power Mobile Information Systems 5 consumption is more suitable for measuring channel quality in a mmWave link, which is modeled by In (16), e i,m (t) is the energy efficiency of the mmWave link i ⟶ m at each interval t, and P RF represents the power consumption coming from an RF chain.
When all the mmWave links determined by the BFT mechanism at each interval t simultaneously employ the common mmWave channel, co-channel interference may occur. For simplicity without loss of generality, we consider the case that the mmWave links (e.g., i ⟶ m and j ⟶ n) are used simultaneously, where the interference power received at GS m at each interval t is estimated by the following formula: In (17), g s j,n ⟶ i,m and g r j,n ⟶ i,m represent the directional transmission gain and directional reception gain between SU j and GS m , respectively. According to formula (1), the directional transmission-reception gain of each path is derived by the following formula: Let ω s j,n ⟶ i,m and ω r j,n ⟶ i,m be the beam offset angle from SU j 's (SU j transmits to GS n ) transmitting beam direction to the position of GS m and that from GS m 's (GS m receives from SU i ) receiving beam direction to the position of SU j , respectively.
We assume that there does not exist multiconnectivity capability for each GSTA. erefore, an GSTA can only connect to one UAV at a time in a UAV-assisted disaster relief network. We set x i,m as an integer variable.
. For convenience, we define a binary variable as follows: For each mmWave link (e.g., i ⟶ m), when considering the interference of all the other mmWave links to it, the signal-to-interference plus noise ratio (SINR) at each interval t is denoted by c sinr i,m , which can be derived by en, the data rate (i.e., throughput) at each interval t of each mmWave link (e.g., i ⟶ m) is given by e sum data rate of all the mmWave links at each interval t is estimated by the following formula: e sum power consumption of all the mmWave links at each interval t is estimated by the following formula: e average system energy efficiency at each interval t is estimated by the following formula:

Radio Channel Allocation Purpose and Problem
Statement. In the paper, under the constraint of a certain number of channels, the purpose of channel allocation is that all the mmWave links obtained by the BFT mechanism should be used concurrently as much as possible to maximize sum data rate. To this end, the optimization problem for channel allocation can be modeled by In (25a) and (25b), χ i,m,c (t) is a binary variable and its value is 1 or 0, where the channel c is used in the mmWave ; C1 is the constraint that any mmWave link between any UAV and any GSTA can only reuse one channel at a certain interval; the difference between r i,m,c (t) and r i,m (t) is only that the former specifies the specific channel while the latter implies the common channel and thus the channel tag is omitted. In essence, r i,m,c (t) is also estimated by formula (21).
It is clear that the optimization problem (25a) and (25b) is the category of combinatorial integer programming problems.
e three-dimensional binary variable χ i,m,c (t) determines the channel selection decisions of each mmWave link i ⟶ m at interval t and thus the space complexity is very large when considering all the mmWave links obtained by the BFT mechanism, all the available channels, and all the intervals. In addition, a centralized optimization technology can theoretically search for the global optimal solution if global information is available. But the side effect is the huge time complexity, especially for combinatorial integer programming problems, such as channel assignment of dynamic network topology. e work in [9] addressed a combinatorial integer programming problem by graph theory in a decentralized manner and proposed the dynamic graph coloring method to solve the channel allocation problem, which reduces both the space complexity and the time complexity. In the research scenario in [9], the construction of the interference graph as the premise of the graph coloring scheme is very simple. However, given the complexity of the scenario in this paper, we firstly need to design a new interference graph construction scheme and then propose a more efficient graph coloring scheme to solve the aforementioned problems with tractable complexity.

Interference Graph Construction and Channel Allocation
In order to transform the optimization problem (25a) and (25b) into a dynamic graph coloring problem, we firstly address interference graph construction problem. We use an undirected graph IG(t) � (L(t), w(t)) to denote the topology structure of interference graph in interval t, where L(t) is the set of the mmWave links obtained by the BFT mechanism in interval t while w(t) is the edge set of interference relationships in interval t between the mmWave links.
In the set L(t), the SINR of each member (e.g., i ⟶ m) in interval t must meet this relation c sinr indicates that an edge exists between i ⟶ m and j ⟶ n, and interference relationship is projected into the edge construction. Also, Θ is used to represent the maximum vertex degree of interference graph IG(t).
In the work in [9], if the desired signal ratio of each user to the mutual interference is below a certain threshold, an edge with respect to interference relationship can be built. In fact, if many links share a common channel, the co-channel interference value of each user is hardly obtained in a simple way, especially in the case of directional mmWave links sharing the same channel. erefore, we propose a new interference graph construction scheme in this section, which is described in Algorithm 1. Since the MUAV can obtain the SNR value of each mmWave link (e.g., i ⟶ m) by the BFT mechanism in each interval t, it estimates the corresponding link energy efficiency e i,m according to formula (16).
Algorithm 1 stores all the links' energy efficiency values into Q � q 1 , . . . , q k , . . . , q K in descending order (see line 4), and thus there is a one-to-one match between the elements of Q and those of e i,m |i ⟶ m ∈ L(t) , where K � Card(L(t)), and Card(·) represents the number of elements in a set. In addition, we set X as a matrix with elements x k K k�1 , and thus there is a one-to-one match between the elements of X and those of x i,m |i ⟶ m ∈ L(t) . ereafter, Algorithm 1 initializes all the elements in X to 0 (see line 5). Next, it chooses the mmWave link from L(t) corresponding to q 1 and then estimates the average system energy efficiency at interval t (i.e., E 1 (t)) based on formula (24) and sets x 1 to 1 (see line 6). Finally, it finds out the corresponding mmWave link based on x 1 and adds it to the vertex set of interference graph (see lines 7∼9). After that, the set AL(t) is transformed from an empty set to the set with one vertex.
Algorithm 1 continues to choose the mmWave link from L(t) corresponding to q 2 and estimates the average system energy efficiency at interval t (i.e., E 2 (t)) (see line 10). If E 2 (t) > ρE 1 (t), Algorithm 1 sets x 2 to 1 (see lines 11∼12), which means that the link corresponding to x 2 can share the same channel as that corresponding to x 1 . Otherwise, it finds out the mmWave link corresponding to x 2 and adds it to the set AL(t) (see lines 14∼16), where its interference relationship with the first member of the set AL(t) is recorded in the set Aw(t) (see line 17). Similarly, subsequent every traversed link that cannot share the same channel with the first link will be added to the set AL(t), and the interference relationship between it and the first link will be recorded in the set Aw(t). After all the links in L(t) have been traversed, Algorithm 2 is invoked to allocate new channels to the edges in the interference graph AIG(t) (see line 19).
After selecting a new channel (see line 1), initializing the temporary set Q' to an empty set (see line 2), and setting the count variable K' and the index variable k to 0 (see lines 3∼4), Algorithm 2 copies the members in the set Q that have not been allocated channels and stores them in the temporary set Q' (see lines 5∼9). en, the members in the set Q' are queued in descending order of their link energy efficiency values (see line 10), and the temporary matrix X' is initialized to 0 in preparation for recording the forthcoming channel allocation (see line 11). Next, Algorithm 2 chooses the mmWave link from L(t) corresponding to q 1 of Q', estimates the average system energy efficiency at interval t (i.e., E 1 (t)) based on formula (24) and sets x 1 to c � 2 (see line 12), and finds out the corresponding mmWave link based on x 1 (see lines 13∼14), where the identified link is marked i ⟶ m and will use channel 2 to communicate.
After determining the first member in the temporary set Q' to use channel 2, Algorithm 2 checks one by one whether other members in the temporary set Q' can share channel 2 with it. If any of other members cannot share channel 2 with the first member, it can be removed from the interference graph AIG(t) (see lines 15∼20). Otherwise, its interference relationship with the first member is added to the edge set Aw(t) of the interference graph AIG(t) (see lines 21∼25).
After processing all the members of the temporary set Q', Algorithm 2 updates the matrix X based on the channel allocation recorded in the matrix X' (see lines 26∼28). en, it checks if there is any record information in the matrix X for the links that have not been allocated channels. If there is any link with unallocated channel and the number of allocated channels has not reached the upper limit, it returns to step 2 to continue channel allocation after selecting a new channel (see lines 29∼31). Otherwise, the channel assignment task is completed.
ere is a constraint factor ρ in both Algorithms 1 and 2, where 0 < ρ < 1. e parameter ρ indicates the degree to which the system may tolerate the loss of system energy efficiency due to mutual interference as the number of concurrent links increases, where the smaller value of ρ will mean the more tolerance of mutual interference but may generate higher number of concurrent links. Because we select the possible concurrent links in descending order in terms of link energy efficiency, the selected first link's energy efficiency (e.g., E 1 (t)) must be the highest one. Subsequently, the energy efficiency of a system consisting of multiple concurrent links must not exceed it. ρE 1 (t) can be regarded as a threshold value which the subsequent system energy efficiency (i.e., E k (t)) will be compared with. Only when the system energy efficiency is higher than this threshold value can the newly added link share the same channel with the original links. In conclusion, the larger value of ρ generates the higher threshold value and thus it is less likely that the newly added link shares the same channel with the original links.

Convergence Performance Analysis
e convergence performance of the proposed channel allocation scheme depends on the time overhead of the BFT process, the interference graph construction algorithm, and the iterative channel allocation algorithm. We adopt the BFT process similar to that in [16], and the main steps are as follows: Initialize the set AL(t) to an empty set (3) Initialize the set Aw(t) to an empty set (4) Store e i,m |i ⟶ m ∈ L(t) into Q � q 1 , . . . , q k , . . . , q K in descending order, where q k � e i,m (5) Initialize X to x 1 � 0, . . . , x k � 0, . . . , x K � 0 (6) Let k be equal to 1, compute E 1 (t), and set x 1 to 1 (7) Get x i,m that matches x 1 (8) Get link i ⟶ m according to x i,m (9) Add i ⟶ m to AL(t) Repeat (10) k � k + 1, and compute E k (t) (11) If E k (t) > ρE 1 (t) then (12) Set x k to 1 (13) Else (14) Get x j,n that matches x k (15) Get link j ⟶ n according to x j,n (16) Add j ⟶ n to AL(t) (17) Add Let c be equal to 2 (2) Initialize the set Q' to an empty set (3) Let K' be equal to 0 (4) Let k be equal to 0 Repeat (5) k � k + 1 (6) If x k ∈ X and x k �� 0 then (7) Take q k out of Q and add it to Q' End if Until k �� K (10) Sort the elements of Q' � q 1 , . . . , q k' , . . . , q K' in descending order (11) Initialize X' to x 1 � 0, . . . , x k' � 0, . . . , x K' � 0 (12) Let k`be equal to 1, compute E k' (t), and set x k' to c (17) Set x k' to c (18) Get x i',m' that matches x k' (19) Get link i' ⟶ m' according to x i',m' (20) Delete i' ⟶ m' and the corresponding edge from AIG(t) (21) Else (22) Get x j,n that matches x k' (23) Get link j ⟶ n according to x j,n (24) Add (i ⟶ m, j ⟶ n) to Aw(t) (25) End if Until k' �� K' (26) If there is any element in X' that has a value of c then (27) Update the value of the corresponding element in X as c (28) End if (29) If there is any element in X that has a value of 0 then (30) If c < C then (31) c � c + 1 and go to 2 (32) Else (33) Go

Mobile Information Systems
(4) e MUAV sends the BFT information of every GSTA to SUAVs through "MUAV FB" frames. (5) Finally, one SUAV will be chosen based on the BFT information to send a "FB" frame to each GSTA. e "FB" frame will notify the GSTA the best beams between itself and its nearby SUAVs.
In Step 1 and Step 2, we assume that the number of training slots is set to K τ and the transmission of each frame takes a training slot size time. In the absence of conflict, the time overhead is K τ + 2 in Step 1 and K τ + 4 in Step 2. So, the time overhead of the entire BFT process is 2K τ + 6 training slots.
In case of conflict, we add a training opportunity. If any conflict occurs in Step 1, the MUAV resends SSW frames in all its beams with conflict signal, where the number of training slots is doubled and thus the time overhead is 2K τ + 2. Similarly, if any conflict occurs in Step 2, each SUAV resends SSW frames in all its beams with conflict signal, where the number of training slots is doubled and thus the time overhead is 2K τ + 4. So, the time overhead of the entire BFT process is 4K τ + 6 training slots. e time overhead of the interference graph construction algorithm is O(K) according to Algorithm 1, while that of the iterative channel allocation algorithm is O(K 2 ) according to Algorithm 2, which is negligible when compared with the BFT time overhead due to the strong processing power of the MUAV.

Experimental Parameter Settings.
In this section, we set up the two simulation scenarios. In the first simulation scenario, as shown in Figure 1, the MUAV hovers at the center of a circular plane with a radius of 100 m, which is 100 m above the ground. Several SUAVs are flying along the edge of the circular plane at the same speed v U . e simulation time is discretized into a series of uniform simulation slices (i.e., a series of uniform intervals). Each simulation time period length is set to an integer multiple of the simulation slices, while each simulation slice is set to an integer multiple of the training slots (i.e., it is more than 4K τ + 6 training slots). e relationship between the three time concepts is shown in Figure 2.
At the beginning of each simulation time period, the speed v U is randomly updated in the range from 0 to 20 km/ h. All the GSATs are randomly distributed in a circular area of the ground with a radius of 150 m centered on the projection point of the MUAV on the ground. e second simulation scenario keeps the whole settings of the first one, but it adds a circular flight trajectory of several SUAVs, which is located 30 meters directly below MUAV with a radius of 50 m. e SUAVs' speed update may be made independently on different flight paths, while the SUAVs' speed update on the same flight path must be synchronized to avoid possible collisions. Figure 3 shows the positional relationship between the two flight paths.
In our simulations, we assume that there is mutual interference caused by multiplexing the same channel if the projections of the UAVs' beam transmission angles on the ground overlap. e main simulation parameters are listed in Table 1.
We use the OMNeT++ network simulation platform (i.e., omnetpp-5.4.1) to implement the above simulation scenario, which runs on the machine with Intel Xeon E5-1680 V4 and Kingston DDR4-2400 16G. During the BFT simulation stage, based on the location information of GSTAs and SUAVs, the simulation module estimates the SNR value of each link according to formula (15). During the interference graph construction and channel allocation simulation stage, the simulation module estimates the energy efficiency value of each link according to formula (16), which will be used as the basis for Algorithms 1 and 2 to determine the order in which all the links are traversed. More details on the simulation results will be analysed below.

Experimental Results and Analysis
. We consider the simulations as shown in  In each simulation experiment, we consider the two scenarios proposed in the above subsection. In each scenario, we further consider three different conditions with the different number of SUAVs. Specifically, for scenario I, the number of SUAVs is set to 2, 4, and 6, respectively; for scenario II, the number of SUAVs is set to 4 (i.e. When the number of GSTAs in the simulation area is fixed at 700 and the parameter K τ is set to 60, Figure 4 shows a comparison of the number of detected GSTAs under the number of training opportunities. From Figure 4, we can see that the average number of detected GSTAs of six scenarios under two training opportunities is about 100 more than that under one training opportunity. However, there is little difference between the three training opportunities and the two training opportunities. erefore, setting two training opportunities can basically find most GSTAs that need communication services and thus adopted in the following simulation experiments. We fix the number of GSTAs in the simulation area to 700 and the number of available channels to 2 and simulate the performance variation trend of the proposed scheme with the number of training slots under six scenarios. e simulation results are shown in Figures 5 to 7.
As can be seen from Figure 5, with the increase of the number of training slots, the convergence time of the proposed scheme in the six scenarios first gradually increases and then shows an increasing trend after a sudden decrease. e main reasons are as follows. On the one hand, the larger number of training slots naturally leads to the longer convergence time; on the other hand, when the number of training slots is very small, the conflict is difficult to avoid, so it is necessary to continue the second training and thus extend the convergence time. However, when the number of training slots is increased beyond a certain value, conflicts no longer occur, and thus the second training is not needed. So the convergence time has a sudden drop, but it increases with the increase of the number of training slots subsequently.
At the same time, we also observe from Figure 5 that there are differences in the performance of the proposed scheme in different scenarios within a certain range of training slots. e reason for this phenomenon is that the number of GSTAs in each UAV's core coverage area is different in different scenarios. For example, for the same number of GSTAs, in 8 SUAVs in scenario II, the average number of the GSTAs competing for each UAV's training slot is the smallest, so the probability of conflict is the lowest. Figure 6 shows that, with the increase of the number of training slots, the average data rate per mmWave link of the proposed scheme in the six scenarios first gradually decreases and then shows a decreasing trend after a sudden increase.
e reason for this phenomenon is somewhat related to the reason for the phenomenon in Figure 5. When the interval of BFT operations is fixed, the longer average convergence time will lead to the smaller data transmission time, and thus less data is received in each simulation slice (i.e., less average data rate per mmWave link). erefore, Figure 6 shows an almost opposite trend to Figure 5. Figure 7 shows the average system energy efficiency of the proposed scheme in the six scenarios under the number of training slots. As the number of training slots increases, the average system energy efficiency also increases. is is because as the number of training slots increases, the more GSTAs users are more likely to be spotted by the UAVs and thus more likely to pick out the concurrent links with higher energy efficiency. However, it is not significant for the increment of average system energy efficiency. Hence, the increasing trend shown in Figure 7 is not obvious.
At the same time, we can also see from Figure 7 that the average energy efficiency per mmWave link is also better in the scenario with the denser UAVs.
is is because the number of concurrent links with the shorter distance and the better position is larger.
We fix the number of training slots to 60 and the number of available channels to 2 and simulate the performance variation trend of the proposed scheme with the number of GSTAs in the simulation area under six scenarios. e simulation results are shown in Figures 8 to 11. Figure 8 shows that the number of detected GSTAs of the proposed scheme in the six scenarios increases with the increase of the number of GSTAs in the simulation area. is means that each UAV has sufficient capability to detect GSTAs when the number of training slots is set to 60.
However, in a scenario with the smaller number of UAVs (e.g., 2 SUAVs in scenario I), the number of detected GSTAs increases less and less as the number of actual GSTAs increases. is is because the actual number of GSTAs in the area covered by these two UAVs is limited. e reason for the faster increase in other scenarios is that there are more UAVs and the discovery capability of each UAV (i.e., enough training slots).
In addition, we also found that, in 6 SUAVs in scenario I, the largest number of GSTAs could be found. is is because the number of UAVs with higher altitude is the largest in this scenario, which enables more GSTAs to be covered with the same beam width. Figure 9 shows that there is no difference in the convergence time for each scenario when the number of GSTAs is very small. is is because the conflict does not occur and the convergence time only depends on a fixed number of training slots. After the number of GSTAs exceeds a certain threshold, the convergence time in each scenario is firstly increased by different amplitude and then stabilized to a specific value. e main reason is that there are conflicts in each scenario, leading to the second training and thus extending convergence time. However, at first, in the different scenarios, not every UAV in every scenario is involved in conflict, so there is a difference in convergence time between them. Finally, with the further increase of the number of GSTAs to 900, each UAV in each scenario has conflicts, which can be solved by adding one training process. erefore, their convergence time tends to be the same. Figure 10 shows that the average data rate per mmWave link is stable at a high level when the number of GSTAs is very small, but presents a sudden downward trend with the increase of the number of GSTAs. As mentioned above, when the interval of BFT operations is fixed, the longer average convergence time will lead to the smaller data transmission time, and thus less data is received in each simulation slice. Based on this, from the trend in Figure 9, we can easily find the reason. When the number of GSTAs is very small, the convergence time is relatively short, and thus the data transmission time is relatively long. When the number of GSTAs exceeds a certain threshold, there is a    sudden increase in the convergence time, so the average data rate per mmWave link presents a sudden downward trend.
Similar to the trend in Figures 7 and 11, it shows that the average energy efficiency per mmWave link of the proposed scheme in the six scenarios increases with the increase of the number of GSTAs in the simulation area. Just as the higher number of training slots can increase the number of the detected GSTAs, the higher number of GSTAs in the simulation area can improve the ability of each UAV's identifying GSTAs. erefore, Figures 7 and 11 and have a similar trend.
We fix the number of training slots to 60 and the number of GSTAs in the simulation area to 700 and simulate the performance variation trend of the proposed scheme with the number of available channels under six scenarios. e simulation results are shown in Figures 12 to 15 . Figure 12 shows that the number of concurrent mmWave links firstly increases with the increase of the number of available channels and then remains at the maximum level when the number of available channels is large enough. is is because the number of RF chains of each UAV is fixed to 8, which means each UAV can communicate well with up to 8 terminals on the ground simultaneously. When the number of available channels is very small, the higher number of available channels is, the      higher number of GSTAs can be served simultaneously by one UAV. When the number of available channels exceeds a certain threshold, each UAV's communication capacity is up to the limit, so the number of concurrent mmWave links remains stable. Figure 13 shows a positive relationship between the number of available channels and the sum data rate. It makes sense because as the number of available channels increases, the number of concurrent mmWave links also increases, which absolutely has a positive impact on the sum data rate for the whole communication system. erefore, Figure 13 shows a similar trend to Figure 12. Figure 14 shows that the average data rate per mmWave link of the proposed scheme in the six scenarios is almost stable with the increase of the number of available channels. Specifically, in 2 SUVA in scenario I, 4 SUAV in scenario I, or 6 SUAV in scenario I, there is no significant difference for the average data rate per mmWave link, no matter what the number of available channels is. at is because the mutual interference between different mmWave links with different UAVs is limited to an extremely low threshold value via the proposed channel allocation approach.
us, the average data rate per mmWave link can remain at a high level. is phenomenon and explanation also apply to the three schemes in scenario II. Figure 15 shows a negative relationship between the number of available channels and the average system energy efficiency of the proposed scheme in the six scenarios. When the number of available channels is very small, the average system energy efficiency is relatively high. When the number of available channels increases, the average system energy efficiency presents a gradual downward trend. As mentioned    above, when the number of available channels increases, the number of served GSTAs (i.e., the number of concurrent mmWave links) also increases; and then the communication interferences in the whole system would increase to a certain level, causing more power consumption. When the number of available channels is large enough, the number of concurrent mmWave links is stable due to the limited number of RF chains of each UAV. erefore, the average system energy efficiency tends to be the same.
Finally, we keep the other parameters unchanged and fix the number of available channels to 12 to show how the variation of the constraint factor ρ can influence the average data rate per mmWave link in our proposed scheme. e simulation result is shown in Figure 16. It can be observed from Figure 16 that the average date rate per mmWave link increases with the increment of the value of ρ in all the six scenarios. e reason is that the tolerance of the system for energy efficiency loss decreases with the increase of ρ, which means that fewer mmWave links can reuse the same channel. Hence, the less co-channel interference results in the higher average data rate. In addition, the widely accepted definition of data rate is the number of successfully transmitted bits per unit of time, which is the opposite of data loss rate. The number of available channels

Conclusions
In this paper, we have investigated the channel allocation problem in integrated mmWave/sub-6 GHz UAV-assisted disaster relief networks. By considering the characteristics of mmWave, we employed the BFT mechanism to find all possible links between ground terminals and UAVs. en, based on the result of BFT process, an interference graph construction algorithm (i.e., Algorithm 1) is proposed to identify the links that cannot be used simultaneously in the same channel. Finally, according to the result of Algorithm 1, an iterative channel allocation algorithm (i.e., Algorithm 2) was designed to allocate the reasonable number of channels for as many communication links as possible with the mutual interference being limited at a threshold value. e simulation results demonstrate that the proposed scheme can achieve an excellent performance in the aspect of average convergence time, average data rate per mmWave link, and average system energy efficiency.
However, there are also some limitations for our work. Firstly, the setting of the three durations (i.e., training slot, simulation slice, and simulation time period) is actually not optimal. Take the duration per simulation slice as an example. If the duration is too short, the BFT process will be executed frequently even though the position of UAVs does not change much, which will lead to unnecessary BFT running overhead. However, if the duration is too long, the BFTprocess will be executed less frequently, which will cause some mmWave links to occupy the allocated channels for a longer time even though the GSTAs are out of the UAVs' coverage area. erefore, the proper BFT running interval should be adaptive to environmental changes. Secondly, the  proposed channel allocation method highly depends on the result of BFT process. If something goes wrong with the BFT mechanism, the channel allocation process will also go wrong. Finally, how to recharge the UAVs is not considered in this paper. How would the UAVs be programmed to recharge? How to deploy recharging stations? What is the impact of them on the model in this paper? ese issues will be considered in our future research plans. In addition, we will also explore the influence of Doppler effect on our model if the high-speed moving scenario is considered.

Data Availability
e simulation data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this study.