Ground-to-Air Communications Beyond 5G: Coordinated Multi-Point Transmission Based on Poisson-Delaunay Triangulation

Based on system-level simulations or measurement trials, the feasibility of communication connectivity for UAVs by LTE networks was evaluated in [4]. Considering 5G technology, the potential use of massive multiple-input multiple-output (MIMO) antennas to support cellular-connected UAVs co-existing with terrestrial UEs was studied in [5]. The theory of stochastic geometry was widely exploited to characterize the performance of cellular-connected UAVs. For instance, the downlink coverage probability in a network including static UAVs and terrestrial UEs was analyzed in [6]. The downlink coverage probability for UAVs in vertical heterogeneous networks was analyzed in [7], where the coverage probability was disclosed as highly dependent on the UAV’s height. Assuming line-of-sight links, an exact coverage probability of cellular-connected static UAVs was derived in [8]. It also reveals that noise and non-line-of-sight links have a negligible effect on the performance of UAVs compared with terrestrial UEs.

On the other hand, the optimization theory was widely applied to acquire optimal performance of cellular-connected UAVs. For instance, the energy efficiency of the rate-splitting multiple access and non-orthogonal multiple access (NOMA) schemes for UAVs was investigated in [9], and the precoding vectors of both methods were optimized accordingly. In [10], the weighted sum rate for both UAVs and terrestrial UEs was maximized by jointly optimizing the UAV’s uplink cell associations and the power allocations. Recently, concerning NOMA, the uplink precoding for UAVs was optimized in [11]. Furthermore, to effectively mitigate inter-cell interference (ICI) to improve the connectivity of UAVs, coordinated multi-point (CoMP) transmission technology comes into sight. In particular, a new cooperative interference cancellation strategy for UAV uplink transmission was proposed in [12]. A closed-form upper bound on the coverage probability of a UAV was derived in [13]. Most recently, a CoMP scheme based on binomial-Delaunay tetrahedralization was proposed in [14], where four aerial BSs serve each UAV within the corresponding tetrahedral cells simultaneously. However, these papers mentioned above do not account for the mobility of UAVs for ease of mathematical tractability.

Considering UAVs’ mobility, a trajectory scheme was designed in [15], subject to the connectivity constraints, including the minimum received signal-to-noise ratio (SNR), the maximum mobility speed, and the initial and final locations. An iterative sub-channel allocation and speed optimization algorithm was developed to maximize the uplink sum rate in [16]. To minimize its mission completion time, the optimal trajectory design for UAVs was studied in [17]. Moreover, to reduce the communication latency of UAVs and their interference with terrestrial UEs, an interference-aware path planning scheme was designed in [18], using machine learning. To minimize the total outage duration, a new trajectory design based on deep reinforcement learning was performed in [19]. Recently, to rapidly model and evaluate the performance of mobile UAVs, two classical spatial mixed mobile models for UAVs with or without pause time were proposed in [20] and [21], respectively, inspired by the classical plane mobile model [22]. Moreover, using the theory of stochastic geometry, the coverage probability for terrestrial UEs and UAVs was analytically derived in [20] and [21], respectively. In this paper, to jointly consider the tractability of the network model and the mobility of UAVs, the classical model, widely known as the modified random waypoint model, will be adopted.

Once the mobility of UAVs is considered, their handoff performance is imperative for the designed network. Regarding two-dimensional (2D) planar scenarios, classical handoff algorithms are usually based on the maximum received power of users [23] or on the user’s dwell-time threshold [24]. There are also speed-based handoff algorithms that can estimate the user’s speed and direction to determine whether a handoff occurs or not [25]. Other handoff criteria include signal power strength and mobile user speed, see, e.g., [26]. For ease of mathematical tractability for handoff performance evaluation, some approximation methods were introduced to simplify optimization problems to get optimal solutions in closed form. For example, given the sojourn time distribution, an effective mobile model was proposed, and the user transition probability from different cells was derived in [27]. Meanwhile, given the location of terrestrial BSs, the movement trajectory of UEs, and the channel state information, an upper bound on the system handoff performance was studied in [28].

Regarding three-dimensional (3D) spatial scenarios, there is little literature focusing on UAVs’ handoff performance. Using the Poisson point process (PPP) to describe the distribution of terrestrial BSs, each projected UAV is associated with its nearest BS, and the resulting polygonal boundaries around BSs form a Poisson-Voronoi tessellation [29]. Based on the Poisson-Voronoi tessellation, the mobile handoff probability for UAVs was studied in [22]. Using CoMP, the cooperative handoff performance for the projected UAV inside conventional hexagonal grids was analytically investigated in [21]. However, the assumption of such regular cells does not coincide with the practical deployment of BSs. Insofar as a mathematically tractable cellular-connected UAV communication architecture integrating CoMP and handoffs is still in its infancy.

The motivation of this paper comes from the works [30, 31], where the Poisson-Delaunay triangulation was exploited to model the 2D cellular networks for terrestrial UEs. A similar idea is employed in this paper to provide service for aerial UAVs, which benefits by maximizing the use of existing cellular network infrastructure and providing reliable communication connectivity for UAVs, in addition to terrestrial UEs. In particular, in this paper, the PPP is applied to model the locations of terrestrial BSs, and the Delaunay tessellation is exploited for cellular network modeling, yielding a network model called Poisson-Delaunay triangulation. Moreover, a frequency-allocation scheme is designed to eliminate ICI by using the theory of circle packing, and its performance is evaluated by the high-dimensional probability theory. Compared with the conventional CoMP strategy with three nearest cooperating BSs [32], simulation and numerical results corroborate that the UAV handoff probability in the proposed scheme reduces to about two-thirds of that in [32]. The main reason behind the reduction is that the nearest distance criterion adopted by the former is more sensitive to UAV mobility.

In this paper, we first develop a mathematically tractable cellular-connected UAV communication architecture based on CoMP transmission. Then, the handoff and coverage probabilities for mobile UAVs are analytically derived. Afterward, a practical frequency allocation scheme is designed to eliminate ICI effectively. In summary, the main contributions of this paper are threefold:

1. 1)

   Network modeling: A cellular-connected UAV communication architecture based on Poisson-Delaunay triangulation is proposed, where CoMP transmission for UAVs is applied to enhance communication connectivity. In particular, each UAV is concurrently served by three terrestrial BSs forming a triangular cell.

2. 2)

   Handoff management and coverage probability analysis: Taking advantage of equivalent circle approximation, a simple but effective handoff management mechanism is first developed. Then, the handoff probability and coverage probability with handoffs for mobile UAVs are derived. Numerical and simulation results show that the handoff probability of the proposed scheme is much lower than that of the conventional CoMP scheme with three nearest cooperating BSs.

3. 3)

   Frequency planning strategy: To effectively eliminate ICI, a dynamic frequency allocation algorithm is designed by means of circle packing, and its performance is evaluated using the theory of high-dimensional probability.

To detail the contributions mentioned above, the remainder of this paper is organized as follows. Section II describes the system model. Section III analyzes the performance of a typical UAV in terms of handoff probability and coverage probability with handoffs. Section IV develops a dynamic frequency allocation strategy. Simulation results are presented and discussed in Section V and, finally, Section VI concludes the paper.

Notation: The operator ${\rm round}(\cdot)$ gives the nearest integer of a real number, and $\lfloor\cdot\rfloor$ returns the maximum integer not larger than a real number. The operator $\mathbb{E}(\cdot)$ means mathematical expectation. The notation $\|\bm{x}\|$ and $\bm{x}^{H}$ denote the $\ell_{2}$-norm and Hermitian transpose of vector $\bm{x}$, respectively. The symbol ${n\choose m}\triangleq\frac{n!}{m!\,(n-m)!}$ indicates binomial coefficient, with $n!$ being the factorial of an integer $n$. The Gamma function is defined as $\Gamma(a)\triangleq\int_{0}^{\infty}{x^{a-1}\exp(-x)}\,{\rm d}x$ for $a>0$. The lower incomplete Gamma function is defined as $\gamma(a,x)\triangleq\int_{0}^{x}t^{a-1}\exp(-t){\rm d}t$ and $\Gamma(v,\theta)$ denotes a Gamma distribution with shape parameter $v$ and scale factor $\theta$. The generalized hypergeometric function ${\displaystyle\,{}_{p}F_{q}(a_{1},\cdots,a_{p};b_{1},\cdots,b_{q};x)\triangleq\sum_{n=0}^{\infty}{\frac{(a_{1})_{n}\cdots(a_{p})_{n}}{(b_{1})_{n}\cdots(b_{q})_{n}}}\,{\frac{x^{n}}{n!}}}$, with $(a)_{n}=a(a+1)\cdots(a+n-1)$ if $n>1$ and $(a)_{n}=1$ if $n=0$. Notice that these special functions can be readily computed using built-in functions in regular numerical software, such as Matlab and Mathematica.

As a companion to the 3D scenario shown in Fig. 1, Fig. 2 shows the BSs and projected UAVs on a 2D plane, denoted by ‘$\circ$’ and ‘$\times$’ marks, respectively. If each projected UAV is associated with its nearest BS, the resulting polygonal boundaries form a Poisson-Voronoi tessellation. The dual Poisson-Delaunay triangulation is illustrated in Fig. 2 as the triangles with solid blue boundaries. For more details on the relationship between Poisson-Voronoi tessellation and its dual Poisson-Delaunay triangulation, please refer to [31] and the references therein.

Given the geometric locations of BSs, the Poisson-Delaunay triangulation in Fig. 2 is uniquely determined and can be efficiently constructed using, e.g., the radial sweep or divide-and-conquer algorithm [33, ch. 4]. Then, for each projected UAV, the CoMP cooperation set can be readily determined. As shown in Fig. 2, for the projected ${\rm UAV}_{1}$, it firstly chooses the nearest BS $A$ and the second nearest BS $B$ as its two associated BSs, thus determining an edge of a cooperative triangular cell. There must be two adjacent triangles that share the same edge and form a quadrilateral $\{A,B,C,D\}$. Among the four BSs at the vertices of the quadrilateral, ${\rm UAV}_{1}$ chooses the two nearest BSs and a third BS between the remaining two opposite BSs that is closer to ${\rm UAV}_{1}$ to form the CoMP cooperating set $\{A,B,C\}$. For a typical UAV, the BSs in $\Phi_{0}$ are called the serving BSs, while all the other BSs are regarded as interfering ones.

We employ the 3D mobility model widely used in the literature, e.g., in [21], which is directly extended from the modified random waypoint model on a 2D plane [34]. As shown in Fig. 3, the movement trace of a UAV can be described by an infinite sequence of tuples: $\{(X_{k-1},X_{k},v)\}$, $\forall X_{k}\in\mathbb{R}^{3}$, where $k$ is the movement epoch and $X_{k}$ is the 3D displacement of the UAV at epoch $k$. During the $k^{\mathrm{th}}$ movement epoch, the UAV is initially from the starting waypoint $X_{k-1}$ at height $H_{k-1}$, then it moves at a constant speed ${v}$ towards the target waypoint $X_{k}$ at height $H_{k}$ with horizontal transition distance $\rho_{k}$ and angle $\phi_{k}$ between the movement direction and the horizontal direction. In general, the horizontal transition length $\rho_{k}$ is characterized by the Rayleigh distribution $f_{\rho_{k}}(\mu)$ with mobility parameter $\mu$. The heights $H_{k-1}$ and $H_{k}$ are chosen uniformly from the finite height interval $[H_{1},H_{2}]$. The selection of waypoints is assumed to be independent and identical for each movement epoch, and there is no pause time at these waypoints [35]. Once the UAV reaches $X_{k}$, it repeats the same procedure to find the next target waypoint $X_{k+1}$ at height $H_{k+1}$, and so on. Let $H$ denote the height of the UAV at a steady state; then its PDF is given by [21]

If $H_{1}=0$, then (II-B) reduces to

which is in accordance with [20, Eq. (3)].

Suppose that all BSs transmit signals with the same power level $P$, then the signal power received at a typical UAV from the three BSs in the CoMP set, $\Phi_{0}=\{A_{0},B_{0},C_{0}\}$, can be computed as $S=\left|\sum_{i\in\Phi_{0}}P^{\frac{1}{2}}d^{-\frac{\alpha}{2}}_{i}\bm{h}^{H}_{i}\bm{\omega}_{i}\right|^{2}$, where $d_{i}$ indicates the Euclidean distance from the $i^{{\rm th}}$ BS to a typical UAV (cf. Fig. 1); $\bm{h}_{i}\in\mathbb{C}^{M\times 1}$ represents the complex-valued channel vector from the $i^{\mathrm{th}}$ serving BS to a typical UAV; $\bm{\omega}_{i}\in\mathbb{C}^{M\times 1}$ is the precoder used at the $i^{{\rm th}}$ BS, and $\alpha>2$ refers to the path-loss exponent. Likewise, the interference power can be computed as $I=\sum_{j\in\Phi\setminus\Phi_{0}}Pd^{-\alpha}_{j,0}|\bm{h}^{H}_{j,0}\bm{\omega}_{j}|^{2}$, where the difference set $\Phi\setminus\Phi_{0}$ represents the point process of interfering BSs; $d_{j,0}$ and $\bm{h}_{j,0}$ stand for the distance and channel vector from the $j^{{\rm th}}$ interfering BS to a typical UAV, respectively.

For general purposes, Nakagami-$m$ fading is assumed throughout the paper to capture a large class of fading environments. Also, full downlink channel state information (CSI) is assumed to be available at BSs. Accordingly, the channel-inverse precoder $\bm{\omega}_{i}$ used by the $i^{{\rm th}}$ BS is given by $\bm{\omega}_{i}={\bm{h}_{i}}/{\|\bm{h}_{i}\|}$. Since $\bm{h}_{i}\in\mathbb{C}^{M\times 1}$ is a random vector with each element $h_{k}$, $k\in\{1,2\cdots,M\}$, being complex-valued Gaussian variable $h_{k}\sim\mathcal{CN}(\sqrt{K},1)$¹¹1The parameter $K$ refers to Ricean factor and $\sqrt{K}$ is exactly the mean of $h_{k}$. The relation between the Ricean factor $K$ and the Nakagami shape factor $m$ is specified in [36, Eq. (2.66)]., $|h_{k}|^{2}$ follows the Gamma distribution with probability density function (PDF) $f(x)=\left({m}/{\Omega}\right)^{m}{x^{m-1}}\exp\left(-{m}/{\Omega}x\right)/{\Gamma(m)}$ [36, Eq. (2.67)] where $m\approx(K+1)^{2}/(2K+1)$ and $\Omega=K+1$. Thus, the channel gain $g_{i}\triangleq|\bm{h}^{H}_{i}\bm{\omega}_{i}|^{2},i\in\{1,2,3\}$, approximately follows the Gamma distribution $\Gamma(m_{1},\Omega_{1})$ with $m_{1}=mM$ and $\Omega_{1}=\Omega/m$, which is the sum of $M$ independent and identically distributed Gamma random variables. As for the interference links, the channel gain $g_{j}\triangleq|\bm{h}^{H}_{j,0}\bm{\omega}_{j}|^{2}$, $j\in\Phi\setminus\Phi_{0}$, also follows Gamma distribution $\Gamma(m_{2},\Omega_{2})$, where $m_{2}$ and $\Omega_{2}$ can be determined by the second-order moment matching method [37, Lemma 7].

According to the discussion in the preceding subsection, the received signal-to-interference ratio (SIR) at a typical UAV can be computed as

Then, similar to [38, 39], given a SIR threshold $\gamma$, the coverage probability with handoffs is defined as

where in (4) ${\rm{Hoff}}$ and ${\rm\overline{Hoff}}$ denote the events that a handoff occurs and no handoff occurs, respectively; the factor $\beta\in[0,1]$ is the probability of connection failure due to handoffs, which reflects the system sensitivity to handoffs. In particular, if $\beta=0$, then the coverage of a UAV is not sensitive to handoffs and becomes equal to the traditional coverage probability; if $\beta=1$, however, every handoff leads to an outage. Clearly, the first term on the right-hand-side of (4) refers to the probability of the joint events that a UAV is in coverage and no handoff occurs, while the second term denotes the probability of the joint events that a UAV is in coverage and a handoff occurs, penalized by the cost of handoff. Finally, in (5) $\mathbb{P}_{\rm{H}}\triangleq\Pr\{\rm Hoff\}$ refers to the handoff probability while $\mathbb{P}_{\rm C}\triangleq\Pr\left\{\eta>\gamma\right\}$ denotes the coverage probability without handoff.

Since a handoff involve possible changes of one, two, or three CoMP BSs within a Delaunay triangle (cf. Fig. 2), it is hard to accurately characterize the coverage area of a virtual cell, namely, the geographic region where a UAV associates with the same CoMP set of BSs before any handoff occurs. As a result, this hinders us from deriving the handoff probability, which is the probability that a UAV moves across adjacent virtual cells during one movement period. To proceed with handoff management, we propose to approximate virtual cells using the dual Voronoi cells.

Now, we explain our handoff management strategy from a geometry perspective. For illustration purposes, Fig. 4a plots four virtual cells with distinct colors, each denoting the geographic region where a UAV associates with the same CoMP set of BSs before any handoff occurs. For a UAV moving within each virtual cell, it is always served by the three BSs at the vertices of the corresponding red triangular cell, and no handoff occurs. Otherwise, a handoff between adjacent cells occurs. On the other hand, Fig. 4b shows the Voronoi tessellation (i.e., the dotted polygons) formed by the black square marks that locate at the centers of the circumcircles of the triangular cells in Fig. 4a. It is observed that the dotted polygons in Fig. 4b can reasonably approximate the corresponding shaded regions with different colors, which represent the virtual cells of our handoff management. Moreover, the area distribution of the dotted polygons in Fig. 4b can be easily computed by that of the Delaunay triangular cells in Fig. 4a by using the theory of stochastic geometry. Consequently, the dotted polygons formed by the nearest-neighbor criterion with respect to the black squares shown in Fig. 4b are employed to analyze the handoff probability, as detailed below.

As illustrated in Fig. 5, define $q_{1}$ and $q_{2}$ in $\mathbb{R}^{2}$ as the horizontal projections of two consecutive waypoints $X_{k-1}(\rho_{k-1},\psi_{k-1},H_{k-1})$ and $X_{k}(\rho_{k},\psi_{k},H_{k})$, respectively. The UAV starting from the projected location $q_{1}$ that is at a distance $r$ from the serving BS $A_{0}$, moves at a speed $v$ to $q_{1}^{\prime}$ at a distance ${v}\cos(\phi_{k})$ in a unit time. Suppose that the UAV does not change its direction during any movement epoch, then $q_{1}^{\prime}$ lies in the segment from $q_{1}$ to $q_{2}$. By using the law of cosines, $q_{1}^{\prime}$ is at distance $R=\left(r^{2}+({v}\cos(\phi_{k}))^{2}+2r{v}\cos(\phi_{k})\cos(\psi)\right)^{1/2}$ from the BS $A_{0}$, where $\psi$ denotes the angle between the horizontal direction and the direction from $q_{1}$ to $q_{1}^{\prime}$ (cf. Fig. 5). As a result, a handoff occurs if there is another BS (e.g., $B_{0}$ in Fig. 5) closer than $R$ to the projected UAV at the new location $q_{1}^{\prime}$. Since a typical UAV moves a distance ${v}\cos(\phi_{k})$ after a unit time, the conditional handoff probability can be formulated in the following lemma.

In light of Lemma 1, by taking the average with respect to $\rho_{k}$, $H_{k-1}$, and $H_{k}$, the handoff probability can be computed as

where $p\triangleq H_{k}-H_{k-1}$ denotes the height difference of a UAV between two consecutive waypoints $X_{k-1}(\rho_{k-1},\psi_{k-1},H_{k-1})$ and $X_{k}(\rho_{k},\psi_{k},H_{k})$, $f_{p}(p)\triangleq{(H_{2}-H_{1})-|p|}/{(H_{2}-H_{1})^{2}}$ with $H_{1}-H_{2}<p<H_{2}-H_{1}$, $f_{\rho_{k}}(\rho_{k})\triangleq 2\pi\mu\rho_{k}\exp(-\pi\mu\rho_{k}^{2})$. Note that $p>0$ corresponds to the rising process of a UAV, whereas $p<0$ to the falling process.

According to the way that each UAV chooses its cooperating set (cf. Fig. 2), the distances between each UAV and its three serving BSs can be approximated to the first three nearest ones $\{d_{1},d_{2},d_{3}\}$. Then, the SIR given by (3) reduces to

By using the theory of Palm distribution, the joint distribution of the horizontal distance parameters $r_{1}$, $r_{2}$, and $r_{3}$ is obtained in [42], expressed as

With the PDF of $H$ shown in (II-B) and the joint PDF of $\{r_{1},r_{2},r_{3}\}$ given by (12), the coverage probability of a typical UAV can be formalized as follows.

Substituting (10) and (13) into (5), an upper bound on the coverage probability with handoffs can be derived, as formalized below.

Like [43], if channel fading is not accounted for, it is equivalent to the case that the shape parameter of Nakagami fading goes to infinity, i.e., $m\rightarrow\infty$. In this case, an accurate expression for the coverage probability of a typical UAV can be explicitly derived. More specifically, in the absence of channel fading, the SIR given by (3) reduces to

By using the relationship between the moment generating function and the characteristic function, the characteristic function of $I_{2}\triangleq\sum_{j\in\Phi\setminus\Phi_{0}}d^{-\alpha}_{j,0}$ can be given by

By recalling [44, Eq. (5)], the coverage probability of a typical UAV can be calculated as

where $\Phi_{\eta_{2}^{-1}}(\omega)$ denotes the characteristic function of $\eta_{2}^{-1}$, computed by

As a result, in the absence of channel fading, the coverage probability with handoffs is explicitly given by

On a 2D Euclidean plane, it is widely known that the highest-density lattice packing of circles is the hexagonal packing arrangement [45, Sec. 1.4], where the centers of the circles are arranged in a hexagonal lattice, and each circle is surrounded by six other circles. In this section, we integrate the method of circle packing into our Delaunay tessellation for frequency planning. The main idea is to use a set of circles to tesselate the coverage area, each covering a cluster of adjacent triangles. The total system bandwidth is normalized to be unity and repeatedly used among different clustering triangular cells. In each cluster of triangular cells, various frequency resources are allocated to different cells. Therefore, inter-cell interference arises only in the triangular cells spanning adjacent circles. On the other hand, as the total system bandwidth is allocated on average to all the cells in a cluster, for a fixed circle area, if the circle covers fewer triangular cells, each cell will have a larger size and be allocated more spectrum resources. This frequency allocation scheme conforms to the actual spectrum demand of engineering applications.

To elaborate on the frequency reuse criterion, the effective interference radius, $\varepsilon$, is defined such that each interfering BS using the same frequency band should be outside this region. The spectral efficiency is chosen as a standard performance metric for the network under study, and a predetermined minimum spectral efficiency $\mathcal{R}_{\rm th}$ in nat/sec/Hz shall be fulfilled, i.e.,

Meanwhile, let $c(o,\varepsilon)$ be a circle with radius $\varepsilon$ centered at the origin $o$, $X\triangleq\left|\sum_{i\in\Phi_{0}}{d_{i}}^{-\frac{\alpha}{2}}\|{h}_{i}\|\right|^{2}$, and $Y_{\varepsilon}\triangleq\sum_{j\in{\Phi}\setminus c(o,\varepsilon)}{d_{j,\,0}^{-\alpha}g_{j}}$. It is straightforward that

which implies that a large $\varepsilon$ is desired from the received SIR viewpoint as $\mathbb{E}\left[Y_{\varepsilon}\right]$ decreases with $\varepsilon$. In contrast, as the said circle $c(o,\varepsilon)$ gets larger, there are more triangular cells in a circle, and each cell gets fewer spectrum resources. To resolve this dilemma, the optimal circular radius, $\varepsilon^{\star}$, is determined by allowing the maximum interference, that is, the inequality (24) takes equality:

where $\mathbb{E}\left[Y_{\varepsilon^{\star}}\right]$ can be computed as [44, Lemma 5]

with $\mathbb{E}[g]=m_{2}\Omega_{2}$ for $g\sim\Gamma(m_{2},\Omega_{2})$ and $\bar{H}=(H_{1}+H_{2})/2$ referring to the average height of UAVs. Meanwhile, $\mathbb{E}[\sum_{i=1}^{3}{d_{i}}^{-\alpha}]$ is evaluated by

with $\{\bm{r}:0<r_{1}\leq r_{2}\leq r_{3}<\infty\}$. Next, inserting (IV-B) and (28) into (26), the optimal circular radius $\varepsilon^{\star}$ for a typical UAV can be computed as

By definition of the circle $c(o,\varepsilon^{\star})$, a cluster of triangular cells is formed within this circle. Each cell in a given cluster uses distinct spectrum resources to eliminate intra-cluster interference. Outside the circle, $c(o,\varepsilon^{\star})$, the triangles with overlapping spectrum resources result in co-channel interference for the typical triangular cell. Therefore, the frequency reuse factor can be defined as

where the numerator and denominator represent the area of the effective interfering circle and the expected area of a typical cell, respectively. In particular, the area of the circle $c(o,\varepsilon^{\star})$ is simply computed by $\mathcal{S}_{c(o,\varepsilon^{\star})}=\pi{\varepsilon^{\star}}^{2}$. Furthermore, according to [40, Table 5.11.1] and [29, Thm. 2.9 and Def. 2.12], we get $\mathbb{E}\left[\mathcal{S}_{\rm cell}\right]=1/(2\lambda)$. As a result, the integer-valued frequency reuse factor in (30) can be explicitly calculated as

Apparently, the frequency reuse factor (equivalently, the number of cells that reserve distinct spectral resources and form a cluster) grows in square law with the effective interference radius, i.e., $\delta\propto{\varepsilon^{\star}}^{2}$.

The basic idea of the proposed frequency allocation strategy is as follows. Given a finite plane of BS intensity $\lambda$, the predetermined spectral efficiency threshold is assumed to be $\mathcal{R}_{\rm th}$, and the average height of UAVs is $\bar{H}$. We first construct the Delaunay triangulation and define the corresponding triangular cells by using, e.g., the radial sweep algorithm or divide-and-conquer algorithm [33, ch. 4]. Then, we compute the radius of circle $\varepsilon^{\star}$ as per (29) and then fill in the target plane with circles of radius $\varepsilon^{\star}$ by using circle packing. Given that there are $k_{1}$ cells in a reference circle say Circle A, the bandwidth reserved for each cell is simply $1/k_{1}$, where the total bandwidth for the whole system is normalized to unity. Finally, evenly assign spectrum resources for each cell according to the number of cells within each circle, starting with Circle A. In each circle, a unique frequency band is assigned to different triangular cells. Meanwhile, the effective interference radius $\varepsilon^{\star}$ is always satisfied for each circle. In summary, the proposed frequency allocation scheme is formalized in Algorithm 1.

For illustration purposes, Fig. 6a shows $7$ circles and $9$ shaded triangular cells sharing the same frequency band. In particular, Circle A with 30 cells is marked as the typical circle, and the red shaded area within Circle A denotes the typical triangular cell with frequency range $[0,1/30]$. Suppose the bandwidth of any cells in interfering circles overlaps with that of the typical cell. In that case, they are called the interfering cells, and the ratio of overlapping bandwidth is taken as a weight value to reflect the strength of interference. For instance, Circle B with 30 cells and Circle C with 33 cells stand for two interfering circles, where the blue shaded cells are the interfering ones. The corresponding overlapping bandwidth is $1/30$ for Circle B, $1/33$ for one cell and $1/30-1/33=1/330$ for another cell in Circle C. Clearly, the sum of the overlapping bandwidth of the two interfering cells in Circle C equals the bandwidth of the typical red cell in Circle A.

For a more intuitive understanding, Fig. 6b gives a schematic diagram of frequency allocation. As said before, assuming the first cell A1 in Circle A to be the typical cell, the cell B1 in Circle B and the cell C1 and part of the cell C2 in Circle C have an overlapping spectrum with the typical cell A1. Therefore, cells B1, C1, and C2 introduce ICI to cell A1. To distinguish the proportions of the overlapping spectrum of different cells in the same Circle, blue and yellow shaded regions are used to mark the corresponding overlapping spectrum regions in Circle C.

On the other hand, Fig. 6a also discloses that there are gaps between tangential filling circles in the proposed circle-filling method. To evaluate the efficiency of the circle packing strategy, we are concerned with the smallest number of circles if seamless coverage is necessary. Suppose that we are given a subset $\mathcal{K}\subset\mathbb{R}^{2}$ and asked to seamlessly cover it by circles of a given radius $\varepsilon>0$. Assume the subset $\mathcal{K}$ can be seamlessly covered by $N_{C}$ circles with radii $\varepsilon$. Comparing their areas, we have $|\mathcal{K}|\leq N_{C}|\mathcal{C}(0,\varepsilon)|$, where $|\cdot|$ denotes the area in $\mathbb{R}^{2}$ and $\mathcal{C}(0,\varepsilon)$ is a Euclidean circle centered at the origin with radius $\varepsilon$. Dividing both sides by $|\mathcal{C}(0,\varepsilon)|$ yields the lower bound $N_{C}\geq|\mathcal{K}|/|\mathcal{C}(0,\varepsilon)|$. On the other hand, we can construct $N_{P}$ closed disjoint circles $\mathcal{C}(x_{i},\varepsilon/2)$ with centers $x_{i}\in\mathcal{K}$ and radii $\varepsilon/2$. While these circles may not need to fit entirely into $\mathcal{K}$, they do fit into a slightly inflated set, namely $\mathcal{K}+\mathcal{C}(0,\varepsilon/2)$. Comparing their areas, we obtain $N_{P}|\mathcal{C}(0,\varepsilon/2)|\leq|\mathcal{K}+\mathcal{C}(0,\varepsilon/2)|$, which leads to the upper bound $N_{P}\leq|\mathcal{K}+\mathcal{C}(0,\varepsilon/2)|/|\mathcal{C}(0,\varepsilon/2)|$. Finally, by recalling [46, Lemma 4.2.8], we have $N_{C}\leq N_{P}$, yielding the following theorem.

Suppose that we fill the unit circle $\mathcal{C}(0,1)$ with different radii ($\varepsilon$) of filling circles. The upper and lower bounds on the number of filling circles for seamless coverage, computed by (33), are listed in Table I, compared with the simulation results of the number of filling circles for the seamless coverage and the number of filling circles for the highest-density lattice packing described in Section IV-A. As expected, the simulation results pertaining to the seamless coverage range between the upper and lower bounds. However, when the radius of filling circles is small (e.g., $\varepsilon=0.025$ or $\varepsilon=0.05$ in Table I), the number of filling circles for the lattice packing is even smaller than the lower bound computed by (33). This observation is of no surprise because the highest-density lattice packing allows a bit of gap among filling circles (rather than seamless coverage), as illustrated in Fig. 6a. When the radius of filling circles becomes larger, the gap becomes smaller such that the numbers of filling circles for the seamless coverage and the lattice packing coincide, as shown in Table I.

After performing the frequency allocation mentioned above, the interfering BSs, not all BSs in $\Phi$, which transmit in the same frequency band, are a thinned version of the original PPP with intensity $\lambda^{\prime}=\lambda/\delta$. Since a thinned version of a PPP is again a PPP, the coverage probability with handoffs can be calculated as per (16) for a typical UAV with $\lambda^{\prime}$ in place of $\lambda$.

Figure 7 depicts the handoff probabilities of UAVs under the proposed Poisson-Delaunay triangulation and the traditional Poisson-Voronoi tessellation with three nearest cooperating BSs, which are illustrated with legends “Delaunay” and “Voronoi: CoMP”, respectively. For comparison purposes, Poisson-Voronoi tessellation without CoMP is also considered, depicted with legend “Voronoi: w/o CoMP”. Finally, the results with legend “Delaunay: Approx.” correspond to the simulation results of handoff probability using the equivalent Voronoi approximation method described in Subsection III-A. Specifically, the left panel of Fig. 7 depicts the handoff probability versus the BS density $\lambda$. It is observed that, for a fixed moving speed of UAVs ($v=20\ {\rm m/s}$), the handoff probability monotonically increases with $\lambda$. Similarly, the right panel of Fig. 7 shows the handoff probability monotonically increases with the moving speed of UAVs for a fixed BS intensity. These observations are because handoffs occur more frequently as the cell size becomes smaller or UAVs fly faster, as expected. On the other hand, the handoff probability of our proposed model is much lower than that of the Voronoi model with three nearest cooperating BSs, especially for large $\lambda$ or high $v$. For instance, it is seen from the right panel of Fig. 7 that, given $\lambda=20\,{\rm BS/km^{2}}$ and $v=40\,{\rm m/s}$, the handoff probabilities of the traditional CoMP scheme and the proposed scheme are about $37\%$ and $24\%$, respectively, thus illustrating the superiority of the proposed scheme. The fundamental reason for the superiority is that, in our scheme, a handoff occurs only if the UAV’s average received power changes, while the distances between a UAV and its three serving BSs are not necessarily the closest. However, in the Poisson-Voronoi network with three nearest BSs, a handoff occurs if any of the three nearest BSs changes, thus resulting in a higher handoff probability.

Figure 7 also shows that the handoff probability of the proposed CoMP scheme is a bit higher than that of the Voronoi scheme without CoMP. This phenomenon is not surprising as the average area of Delaunay triangular cells is half of their dual Voronoi cells [30]. On the other hand, it is observed that the simulation results based on the equivalent Voronoi approximation method match well with the simulation results, demonstrating this method’s feasibility. Finally, it is seen that the numerical results computed as per (10) agree well with the simulation ones, which corroborates the effectiveness of our analysis.

In this subsection, simulation results about the coverage probability with handoffs are firstly depicted, then exact numerical results about the coverage probability in the absence of fading are shown. Next, the effect of varying heights of UAVs is discussed, and the coverage probability after frequency allocation is finally illustrated.

Figure 12 illustrates the coverage probability of a typical UAV versus the SIR threshold $\gamma$ in dB, regarding the proposed Delaunay CoMP scheme and the Poisson-Voronoi tessellation scheme with three nearest cooperating BSs. For comparison purposes, the simulation results of the Poisson-Voronoi tessellation without CoMP are also plotted. If no handoff is accounted for (i.e., $\beta=0$), Fig. 12a shows that the proposed Delaunay scheme achieves a slightly lower coverage probability than the Voronoi scheme with three cooperating BSs. This is because the latter always chooses the three nearest serving BSs through exhaustive searching. However, when handoff is accounted for, for instance, in the case of $\beta=0.5$, Fig. 12b shows that the coverage probability of our proposed Delaunay scheme outperforms that of the Voronoi scheme with three nearest cooperating BSs, especially at low SIR threshold. This comparison concludes that our CoMP scheme is preferable as UAVs inevitably require handoffs.