{'fraction_non_alphanumeric': 0.03920160284277776, 'fraction_numerical': 0.009148300759838204, 'mean_word_length': 4.317820886521258, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Network-Connected UAV: 3D System Modeling and Coverage Performance Analysis

With growing popularity, unmanned aerial vehicles (UAVs) are pivotally extending conventional terrestrial Internet of Things (IoT) into the sky. To enable high-performance two-way communications of UAVs with their ground pilots/users, cellular network-connected UAV has drawn significant interests recently. Among others, an important issue is whether the existing cellular network, designed mainly for terrestrial users, is also able to effectively cover the new UAV users in the three-dimensional (3D) space for both uplink and downlink communications. Such 3D coverage analysis is challenging due to the unique air-ground channel characteristics, the resulted interference issue with terrestrial communication, and the non-uniform 3D antenna gain pattern of ground base station (GBS) in practice. Particularly, high-altitude UAV often possesses a high probability of line-of-sight (LoS) channels with a large number of GBSs, while their random binary (LoS/Non-LoS) channel states and (on/off) activities give rise to exponentially large number of discrete UAV-GBS association/interference states, rendering coverage analysis more difficult. This paper presents a new 3D system model to incorporate UAV users and proposes an analytical framework to characterize their uplink/downlink 3D coverage performance. To tackle the above exponential complexity, we introduce a generalized Poisson multinomial (GPM) distribution to model the discrete interference states, and a novel lattice approximation (LA) technique to approximate the non-lattice GPM variable and obtain the interference distribution efficiently with high accuracy. The 3D coverage analysis is validated by extensive numerical results, which also show effects of key system parameters such as cell loading factor, GBS antenna downtilt, UAV altitude and antenna beamwidth.


I. INTRODUCTION
With enhanced functionality and reducing cost, unmanned aerial vehicles (UAVs), or so-called drones, have found fastgrowing applications over recent years in the civilian domain such as for cargo delivery, precise agriculture, aerial imaging, search and rescue, etc. In particular, UAV can be employed as aerial communication platform [1] to provide wireless connectivity for the ground users and Internet of Things (IoT) devices [2] when traditional terrestrial networks are unavailable, insufficient or costly to deploy. Typical applications include UAVaided communication offloading for temporary hotspot regions [3]- [9]; mobile data relaying between distant ground users [10]; and efficient information dissemination or data collection in IoT and sensor networks [11] [12], etc. On the other hand, UAVs can be integrated into the existing and future cellular networks as new aerial user equipments (UEs) to enable their two-way communications with the terrestrial users efficiently, thus extending the IoT to the sky, known as the Internet of Drones (IoD) [13]. To achieve the cellular-enabled IoD, it is of paramount importance to ensure that all UAVs can operate safely and reliably, even in harsh environments. This calls for ultra-reliable, low-latency, and secure communication links between the ground base stations (GBSs) and the UAV for supporting the critical control and non-payload communications (CNPC). Moreover, in many applications such as video streaming, surveillance and aerial imaging, UAVs generally require high-capacity data communication links with the GBSs so as to timely send the payload data (such as high-quality images and videos) to the end users.
To enable high-performance two-way communications between UAVs and ground users, the existing 4G (fourthgeneration) LTE (Long Term Evolution) or forthcoming 5G (fifth-generation) cellular networks can be leveraged thanks to their almost ubiquitous accessibility and superior performance. As a result, network-connected UAV communications have drawn significant interests recently (see e.g. [14]- [17] and the references therein). In fact, the 3rd Generation Partnership Project (3GPP) had launched a new work item to investigate the various issues and their solutions for UAV communications in the current LTE network [18]. Moreover, increasingly more field trials have been conducted on using terrestrial cellular networks to provide wireless connectivity for UAVs [19] [20].
Among others, one critical issue to address for cellularenabled UAV communications is whether the cellular network is able to provide reliable three-dimensional (3D) coverage for the UAVs at various altitude in both the uplink and downlink communications 1 . To achieve reliable UAV-GBS communications, the uplink or downlink signal-to-noise ratio (SNR), with co-channel interference treated as additional noise, needs to be no smaller than a predefined threshold; otherwise, an outage will occur. For the UAV's uplink, it may transmit multimedia data to the GBS, which requires high data rate and correspondingly high SNR threshold. On the other hand, the UAV's downlink communication needs to support the critical CNPC data from its associated GBS, which is typically of lower data rate but requires higher reliability (lower outage probability) compared to the uplink. The 3D coverage probability of the UAV uplink/downlink communication is thus defined as the corresponding average non-outage probability of a UAV uniformly located in a given 3D space.
The 3D coverage performance analysis for aerial users is considerably different from its two-dimensional (2D) counterpart for ground UEs in the traditional cellular network [21], due to their distinct channel characteristics and resulted interference effects. Specifically, for ground UEs, their channels with GBSs usually exhibit severe pathloss and prominent small-scale fading due to the rich scattering environment especially in urban areas. In contrast, for UAVs in the sky far above the GBSs, their communication signals are more likely to propagate through free space with few obstacles, and hence line-ofsight (LoS) links usually exist with a high probability, which increases with the UAV altitude in general [22]. Although LoS channels entail more reliable communication between the UAV and its serving GBS as compared to the terrestrial UEs, they also cause more severe uplink/downlink interference to/from a larger number of non-associated GBSs. This thus calls for effective interference management techniques such as multi-cell coordinated channel assignment and transmission [23], advanced antenna beamforming techniques at the UAV [24] and GBS [25], interference-aware UAV path planning [17], [26], etc. Besides the LoS dominant air-ground channel, the 3D GBS antenna pattern also has a significant effect on the UAV's coverage performance. The commonly adopted GBS antenna pattern in 3GPP [27] or existing literature [21] [28] is usually simplified in the vertical domain, e.g., by specifying only fixed gains for the antenna mainlobe and sidelobes, respectively. In practice, however, the GBS antenna is usually tilted downward to support ground UEs [19], and hence likely to communicate with aerial UEs in its sidelobes only. Therefore, a more practically refined model for the GBS antenna pattern is needed to characterize the 3D coverage performance of the UAV accurately, especially for the case with non-uniform sidelobe gains and even nulls between them.
This paper thus focuses on modeling the cellular-enabled UAV communication system and analyzing its 3D coverage performance for both the uplink and downlink, by taking into account the unique air-ground channel characteristics and the practical non-uniform GBS antenna pattern. Specifically, we consider a practical UAV-GBS association strategy where the UAV is associated with the GBS that provides the strongest channel gain with it, which can be implemented by comparing the reference signal received power (RSRP) of the downlink signals sent by the GBSs. To capture the essential feature of air-ground channel, we adopt a simplified but practical binary channel state model, comprising only the two states of LoS or non-LoS (NLoS) with different probabilities of occurrence [22]. As a result, the UAV might be associated with a distant instead of nearby GBS, depending on its channel states with the GBSs as well as angles with their antennas. Besides GBS association, the coverage performance also depends on the interference with co-channel GBSs, which is a discrete random variable (RV) (as opposed to continuous RV in terrestrial communication) due to the probabilistic LoS/NLoS channel model and the random on/off activities of co-channel GBSs. As a result, the 3D coverage analysis invokes discrete channel and interference states, and their numbers increase exponentially with the number of involved GBSs within the UAV's signal coverage, which is practically large due to the high probability of LoS. To our best knowledge, the coverage performance of the cellular network under such large discrete channel/interference states has not been addressed, which is fundamentally different from that of the terrestrial 2D network with continuous fading channel/interference [21], or the 3D air-ground network [6] [9] without considering the probabilistic LoS/NLoS channel. The main contributions of this paper are summarized as follows.
• First, we present a 3D system model for the UAV-GBS uplink/downlink communications, which includes the cellular network model, the 3D patterns of the GBS and UAV antennas, and the 3D air-ground channel. Note that our model is applicable to any given 3D patterns of the GBS/UAV antennas. • Second, we propose an analytical framework to characterize the uplink/downlink 3D coverage (average nonoutage probability) performance of the UAVs. The new contributions mainly include the consideration of the probabilistic LoS/NLoS channel states, the UAV-GBS association and coupled downlink interference analysis, and the resulted discrete SNR distribution characterization. To this end, an efficient sorting algorithm is proposed to analyze the UAV-GBS association and uplink SNR distribution, which significantly reduces the complexity from an exponential order with the number of involved GBSs by exhaustive search to a linear order. Moreover, to resolve the coupling between the UAV-GBS association and downlink interference, the downlink SNR distribution characterization is reduced to deriving the interference distribution conditioned on a given UAV-GBS association. • Third, we model the conditional downlink discrete interference given the associated GBS as a new distribution termed as generalized Poisson multinomial (GPM), which take into account the LoS/NLoS channel states and random on/off channel activities of all co-channel GBSs. However, the size of the sample space of the discrete GPM RV increases exponentially with the number of cochannel GBSs, and as a result its cumulative distribution function (cdf) is prohibitive to compute via the bruteforce enumeration-based method. Furthermore, the conditional interference distribution needs to be evaluated over all possible associated GBSs to obtain the downlink SNR distribution and hence the coverage probability. To reduce such high complexity, we propose a new and efficient method to obtain the conditional interference distribution with high accuracy, named the lattice approximation (LA) method, which converts a general non-lattice distributed GPM RV into a lattice distribution with a bounded size of the value space, and then applies the efficient fast Fourier transform (FFT) on its characteristic function to obtain the approximate interference distribution with high accuracy. • Finally, extensive numerical results are provided to validate our analysis and reveal insights for system design. First, it is shown that the GBS antenna pattern has a significant impact on the spatial distribution of non-outage probability, which is useful for UAV path planning/movement control. Moreover, a large downtilt angle leads to overall smaller GBS antenna gain for the UAV above the GBS height, which affects the uplink and downlink coverage probabilities in different ways.
Second, as the UAV altitude increases, the UAV-GBS link distance increases while so does the LoS probability, both of which affect the link strength and hence the coverage probability, but in opposite ways. In addition, a high LoS probability at high altitude leads to small variation of the coverage probability versus the SNR threshold. Third, a high network loading factor with more terrestrial UEs reduces the UAV coverage probability as well as its variation. Finally, besides the antenna gain, applying directional antenna at the UAV limits the coverage range of its antenna mainlobe and hence the number of covered GBSs, which effectively mitigates the interference at high UAV altitude and thus improves the coverage performance. The rest of this paper is organized as follows. Section II introduces the 3D system model. The 3D coverage analysis for the uplink/downlink communication is presented in Section III and Section IV, respectively. Numerical results are provided in Section V, followed by conclusions in Section VI.
Notations: R denotes the set of real numbers; Z denotes the set of integer numbers; P{·} denotes the probability of an event; E{·} denotes the expectation of an RV; · denotes the Euclidean norm; | · | denotes the cardinality of a set; \· denotes the set minus operation; denotes the set union; denotes the set intersection; and ∅ denotes the empty set.

II. SYSTEM MODEL A. Cellular Network Model
Consider a cellular network with the classic hexagonal grid cell layout 2 where each GBS is at the center of its cell with the inter-cell distance of D meters (m) and a given frequency reuse factor of ρ = 1/F, where F ≥ 3 and F ∈ Z. An example of the network with ρ = 1/3 is shown in Fig. 1(a), where the set of GBSs are denoted as B = {0, 1, 2, · · · } and represented by circles of different colors. The whole spectrum is equally divided into F = 3 orthogonal bands, each of which is reused by the GBSs of the same color, denoted by the co-channel GBS sets B 1 (yellow), B 2 (blue) and B 3 (green), respectively, e.g., B 1 = {1, 3, 5, 9, 13, 17, · · · }. Note that generally B f , f = 1, · · · , F are orthogonal and F f=1 B f = B. For simplicity, we assume that all GBSs are at an identical height of H b m, and the 3D coordinate of GBS i ∈ B is denoted by Consider a typical UAV user flying at an altitude of H u m, with H u > H b , whose 3D coordinate is denoted by u (x u , y u , H u ), where x u , y u ∈ R. Without loss of generality, we assume that GBS 0 is at the origin and the UAV's horizontal position is randomly located inside the reference cell 0.
In the conventional cellular network, each ground UE is usually associated with one of its nearby GBSs for communication based on the average channel gain mainly determined by the distance-dependent pathloss and shadowing. However, for the UAV user flying at high altitude, it is possible that the UAV connects to a distant GBS, due to the random realization of LoS/NLoS links, as well as the variation of the GBS 2 Note that although we use the regular grid as an example, our analysis can be extended to arbitrary cell topology in practice.  antenna radiation pattern in the elevation domain, especially the sidelobes and the nulls in between. An illustrative example is given in Fig. 2, where the UAV user is served by the more distant GBS 3 via a strong antenna sidelobe instead of the nearby GBS 1 with a possible null of antenna gain. Denote i s ∈ B as the index of the GBS that serves the UAV of interest, and f s as the index of the set of co-channel GBSs of GBS i s , i.e., i s ∈ B fs . For example, suppose that the serving GBS is i s = 3, then the set of co-channel GBSs is B fs = B 1 as represented by yellow circles in Fig. 1 We consider that orthogonal time-frequency resource blocks (RBs) are assigned to the UAV for its uplink and downlink communications, respectively, by its associated serving GBS i s . We assume that the RBs are assigned to users by the GBSs in each co-channel GBS set B f , independently. For each cochannel GBS i ∈ B fs \ i s , we assume that the RBs assigned to the UAV by GBS i s in the uplink/downlink are simultaneously used to serve a ground UE with probability ω ul,i with 0 < ω ul,i < 1, and ω dl,i with 0 < ω dl,i < 1, respectively. The channel active probabilities ω ul,i and ω dl,i reflect the current uplink/downlink loading factor of GBS i, respectively, which are assumed to be given.
Specifically, we define a binary variable µ i = 1 to indicate that the co-channel GBS i ∈ B fs \i s is active in communication with a ground UE in the same assigned RB as that of the UAV in the downlink, and otherwise µ i = 0. Then the downlink µ i 's for different co-channel GBSs i ∈ B fs \i s are independent Bernoulli RVs with parameter ω dl,i , respectively. Similarly, we can define another binary variable ν i to represent the channel activity of co-channel GBS i ∈ B fs \ i s on the uplink RB used by the UAV.

B. Antenna Model
We consider a practical GBS antenna gain pattern which is omnidirectional in the horizontal plane but vertically directional. 3 Denote θ ∈ (−90 • , 90 • ] as the elevation angle upward from the horizontal plane of the GBS antenna, as shown in Fig. 2. Denote the GBS antenna gain at an elevation angle θ as G b (θ). For a GBS i ∈ B, the elevation angle θ i as seen by the UAV is given by The specific GBS antenna gain pattern depends on the GBS antenna type and configuration. In existing cellular networks, the GBS antenna is usually tilted downward to support ground UEs, where the antenna boresight direction is electrically or mechanically downtilted with an elevation angle θ tilt degree (θ tilt < 0), as shown in Fig. 2. For the purpose of exposition, we consider in this paper the GBS antenna pattern synthesized by a uniform linear array (ULA) [29] with K co-polarized dipole antenna elements placed vertically with d e spacing between elements and electrically steered with downtilt angle θ tilt . According to [29], the power gain pattern of the ULA is given by G e,max cos 2 θ is the power gain pattern of each dipole antenna element with G e,max denoting its maximum value; and J(θ) is the normalized array factor of the ULA with ϑ 2π λ d e (sin θ − sin θ tilt ) in radian (rad) and λ being the wavelength. Note that in the downtilt direction, i.e., θ = θ tilt , we have ϑ = 0 and J(θ tilt ) = √ K, and hence which approximately achieves the maximum antenna gain in (2) for small downtilt angle in practice, e.g., |θ tilt | ≤ 20 • . For illustration, the power gain pattern of the ULA in the 3D space for θ tilt = −10 • is plotted in Fig. 3 (a), where the half-wave dipole element is used with G e,max = 1.64 and the following parameters: K = 10, d e = 0.5λ and carrier frequency f c = 2 GHz. The corresponding 2D patterns in the elevation domain for θ tilt = −10 • and θ tilt = −20 • are plotted in Fig. 3 (b) and (c), respectively. Similar examples of GBS antenna pattern synthesized by antenna arrays can be found in the latest 3GPP technical report [30]. Note that our proposed analytical framework in this paper is general and can be applied to any given GBS antenna pattern.
On the other hand, the power gain pattern of the UAV antenna also plays an important role in the link performance. In particular, the use of directional antenna with different boresight directions at different UAV altitudes can effectively confine the interfered/interfering co-channel GBSs in the uplink/downlink within a limited region and hence alleviate 3 Our analysis can be readily extended to the case with sectorized antenna pattern of the GBS in the horizontal plane. the severe air-ground interference issue due to strong LoS channels. Note that depending on the hardware configuration and flying altitude of the UAV, the boresight direction of the UAV antenna in general could be adjustable within a certain range. For simplicity, in this paper, we consider a given boresight direction of the UAV antenna at each altitude. For the purpose of exposition, we assume that the UAV is equipped with a directional antenna whose boresight direction is pointing downward to the ground, and the azimuth and elevation half-power beamwidths are both 2Φ u degrees (deg) with Φ u ∈ (0 • , 90 • ), as shown in Fig. 2. Furthermore, the corresponding antenna power gain in direction (φ a , φ e ) can be practically approximated as where G 0 = 30000 2 2 = 7500; φ a and φ e denote the azimuth and elevation angles from the antenna boresight direction, respectively [29]. Note that in practice, g 0 satisfies 0 < g 0 G 0 /Φ 2 u , and for simplicity we assume g 0 = 0 in this paper. Further note that, in the special case with Φ u = 90 • , we have G u (φ a , φ e ) ≈ 1, ∀φ a , φ e , which reduces to the case with an isotropic antenna. The antenna gain of the UAV as seen by a GBS i ∈ B is thus given by where is the radius of the coverage area of the UAV antenna mainlobe projected on the horizontal plane at the GBS height, as shown in Fig. 2. As a result, the UAV can only communicate/interfere with the GBSs within its antenna coverage area.

C. Channel Model
Due to the high altitude of the UAV, LoS channel exists with a high probability for the UAV-GBS link, which increases with the UAV altitude in practice [22]. According to [31], the received signals in UAV-ground communications mainly constitute three components, namely LoS signal, strong reflected NLoS signals and multiple reflected signals which cause multi-path fading, each with a certain probability of occurrence. Typically, the probability of receiving LoS and strong NLoS signals are significantly higher than that of multipath fading [32]. Therefore, for simplicity as well as capturing the main characteristic of UAV-ground channels, we ignore the multi-path fading in this paper and consider only the dominant LoS and NLoS components. The authors in [33] proposed a segmented regression approach to further extend the channel model to constitute multiple propagation groups besides LoS/NLoS (e.g., obstructed LoS), which is shown to be able to reconstruct the segmented structure of the UAVground propagation conditions observed in practice and by ray tracing simulations. Note that our proposed analytical framework in this paper can be extended to any finite states of the channel realization.
In practice, the LoS/NLoS probability and associated pathloss depend on the density and height of buildings in a given environment, as well as the relative position between the UAV and GBS. Therefore, in a given environment, the channel power gain between GBS i ∈ B and the UAV at position u can be expressed as where h L,i and p L,i denote the channel power gain and occurrence probability of the LoS channel, respectively, while h NL,i and p NL,i denote the counterparts for the NLoS channel, respectively, and p NL,i = 1 − p L,i . For the specific forms of the involved functions in (6), interested readers may refer to the simplified formula in [5] or the empirical formula in [22]. Furthermore, we define a binary variable δ i to represent the event of LoS/NLoS channel occurrence, i.e., δ i 1, LoS channel between UAV and GBS i; 0, NLoS channel between UAV and GBS i.
Then the δ i 's for different GBSs i ∈ B are independent Bernoulli RVs with parameter p L,i , respectively. Then the channel gain in (6) can be rewritten as

D. UAV-GBS Association
Denote C i (u) as the overall power gain by taking into account both the channel power gain and the antenna power gains of the UAV and GBS i ∈ B, which is given by In this paper, we consider a practical user association rule where the UAV is associated with the GBS that provides the strongest power gain with it, 4 which can be implemented by comparing the RSRP of the downlink beacon signals sent by the GBSs in B. Assume that the beacon signals are sent using the same transmit power. Then the serving GBS i s is selected as the one with the largest power gain, i.e., 4 Here we assume that all GBSs have available RBs to support the uplink/downlink communications of a new UAV user.

E. Uplink Communication
Besides the direct communication link with the associated GBS, the interference links also affect the 3D coverage performance of the UAV. Specifically, there are two new types of interference in a cellular network supporting both ground UEs and UAV users, which are respectively the interference between UAV users and ground UEs, and that among different UAVs even when they are not densely distributed. In this paper, we assume that different UAVs are assigned with orthogonal channels and thus there is no inter-UAV interference. We further assume that the traditional interference issue among ground UEs in co-channel cells is effectively resolved by existing techniques such as cell planning, frequency reuse, dynamic RB allocation, power control, beamforming, etc. Therefore, we focus on characterizing the new air-ground interference between the UAV user and GBSs 5 .
For uplink communication, we assume that the UAV transmits with power P u in Watt (W), which is capped by the maximum transmit power P max . The received SNR 6 at the receiver of the serving GBS i s is thus given by where β 0 P u /σ 2 and the receiver noise is assumed to be additive white Gaussian noise (AWGN) with zero mean and power σ 2 .
Note that the uplink SNR γ ul (u) at each given UAV location u is an RV due to the probabilistic LoS channel occurrence δ i , i ∈ B. Moreover, γ ul (u) also depends on the uplink transmit power P u . We consider that the UAV has an uplink SNR requirement of γ ul (u) ≥ η ul in order to be in non-outage, where η ul > 0 is a pre-defined SNR threshold. The uplink outage probability of the UAV at location u is thus defined as where F γul(u) (·) denotes the cdf 7 of γ ul (u). Assume that the UAV performs power control with the objective to satisfy the SNR requirement, while at the same time reducing the interference to co-channel GBSs, subject to the maximum transmit power P max . Specifically, in the case where the channel gain C is of the serving GBS i s is sufficiently large such that ≥ η ul , then the SNR requirement can be satisfied and the UAV can reduce its transmit power to P u = η ul σ 2 /C is (u). On the other hand, in the case with PmaxCi s (u) σ 2 < η ul , the SNR requirement cannot be satisfied even with the maximum transmit power P max ; as a result, the UAV is said to be in the uplink "coverage hole". In this case, the UAV can keep silent to avoid the interference to other co-channel GBSs. In the rest of this paper, without loss of generality, we consider P u = P max in order to characterize the (maximum) non-outage probability in the uplink.
In practice, if the UAV is allowed to adjust its position u within a certain region, then the UAV can move to other position u with potentially better SNR to get out of the coverage hole. For this purpose, the spatial distribution of the outage probability within the considered region is useful to guide the direction of the UAV movement in practice. In Section V, we will investigate the spatial pattern of outage probability distribution for such applications.
Next, we define the uplink coverage probability of the UAV in the 3D space S ⊂ R 3 as the spatial average of non-outage probability over the space, i.e., where the UAV position u is assumed to be uniformly distributed in the considered 3D space S. More specifically, we consider two subspaces. First, define the 2D subspace at a given UAV altitude H u within the 2D horizontal area A ⊂ R 2 as Second, define the 3D subspace between an altitude range Accordingly, we can investigate the coverage probability at a certain UAV altitude or altitude range in the cellular network.
In the considered hexagonal cell layout, thanks to symmetry, the coverage probability of the UAV can be obtained by applying (13) to the uniformly distributed UAV horizontal positions within the area of the reference cell 0, denoted as A 0 . Note that by symmetry, we can further divide each reference cell into six centrally-symmetric triangular parts, and then we only need to average over one such part A as shown in Fig. 1(a) to compute the coverage probability. Therefore, the overall uplink coverage probability of the UAV at an altitude H u is given by

F. Downlink Communication
For downlink communication, the UAV receives interference from its active co-channel GBSs (see Fig. 1(b)), which cannot be ignored due to the strong LoS channels. Assume that each GBS transmits with the same power P b (W) to the UAV/ground UE in the donwlink. Then the desired signal power received by the UAV from the serving GBS i s is given by P b C is (u); while the interference power from each active co-channel GBS i ∈ B fs \ i s is given by P b I i (u), where I i (u) µ i C i (u). The aggregate interference power (normalized by P b ) received by the UAV is thus given by As a result, the SNR received at the UAV is given by where α 0 σ 2 /P b . The downlink SNR γ dl (u) at UAV location u is an RV due to the random channel activity µ i of co-channel GBSs and the probabilistic LoS channel occurrence We assume that the UAV has a downlink SNR requirement of γ dl (u) ≥ η dl in non-outage, where η dl > 0 is a pre-defined SNR threshold. The downlink outage probability of the UAV at location u can then be defined as where F γdl(u) (·) denotes the cdf of γ dl (u). Based on (20)

III. UAV-GBS ASSOCIATION AND UPLINK OUTAGE ANALYSIS
In this section, we analyze the UAV's uplink outage probability in (12) at any given location u, which is essential to the characterization of the 3D coverage probability in (13). For uplink communication, the SNR γ ul (u) is determined by the UAV-GBS association based on the strongest power gain, which is, however, random due to the probabilistic LoS channel occurrence δ i , i ∈ B.
For a given UAV location u, suppose that the antenna power gains of the UAV and GBS i ∈ B, the channel power gains h L,i (u) and h NL,i (u) of LoS/NLoS channels, and the LoS probability p L,i (u) are all given. For notation simplicity, we drop (u) in the following analysis. Then the uplink SNR in (11) can be rewritten as where each term C i can take two possible values of C L,i G u,i G b θ i h L,i or C NL,i G u,i G b θ i h NL,i , depending on the realization of the Bernoulli RV δ i which takes the value of 1 with probability p L,i or 0 with probability 1 − p L,i . In general, we have C L,i > C NL,i , ∀i ∈ B.
In order to obtain the probability mass function (pmf) 9 of the uplink SNR γ ul , a direct method is to enumerate 2 |B| possible combinations of C i , i ∈ B, then find for each realization the maximum C i and its probability of occurrence. This enumeration-based method has an exponential complexity in terms of |B|, which can be practically large (e.g., when the UAV altitude H u is high and the UAV antenna beamwidth Φ u is large). In order to reduce such complexity, we propose an efficient algorithm to obtain the pmf of γ ul in Algorithm 1. The key idea is to sort the LoS channel gains {C L,i } i∈B in descending order, then obtain the maximum C i and its probability one by one based on the sorted order. Specifically, denote the ordered index of the GBSs as i m , m = 1, · · · , |B|. First, we have P{γ ul = β 0 C L,i1 } = p L,i1 , since max i∈B C i = C L,i1 when the channel realization between the UAV and GBS i 1 is LoS. Similarly, we have P{γ ul = β 0 C L,im } = p L,im m−1 j=1 (1 − p L,ij ), since max i∈B C i = C L,im when the channel realization between the UAV and GBS i m is LoS while those between the UAV and GBSs i 1 , · · · , i m−1 are NLoS. As a result, the worse-case complexity of Algorithm 1 is only linear in the maximum number of iterations |B|.
To further reduce the number of iterations in Algorithm 1, two early stopping criteria can be applied. First, let C max NL max i∈B C NL,i be the maximum NLoS channel gain among the GBSs in B. Then the algorithm can stop early if the current LoS channel gain C L,im is smaller than C max NL , since the maximum C i cannot be smaller than C max NL . Second, in the mth iteration, if the probability term m−1 j=1 (1 − p L,ij ) is lower than a prescribed small threshold value > 0, then we can also neglect the rest of iterations and stop the algorithm early. Finally, after the pmf of uplink SNR γ ul (u) at UAV location u is derived, we can then obtain its cdf and hence the uplink outage probability in (12).

IV. DOWNLINK OUTAGE ANALYSIS
In this section, we analyze the UAV's downlink outage probability in (20) at any given location u, which is more involved as compared to that in the uplink derived in the previous section. For downlink communication, the SNR γ dl (u) depends on not only the UAV-GBS association, but also the aggregate interference I(u) from all co-channel GBSs i ∈ B fs \ i s , which is a discrete RV due to the probabilistic LoS channel occurrence δ i and the random channel activity µ i . First, we derive the general cdf expression of γ dl (u), by resolving the coupling between the UAV-GBS association and the aggregate interference based on conditional probabilities. Then conditioned on the associated GBS, we further investigate the cdf of aggregate interference with practically a large number of co-channel GBSs in our considered system. As a result, the computational complexity of the enumerationbased method to directly compute the discrete conditional interference cdf is exponential in the number of co-channel GBSs, which is prohibitive for implementation. To reduce such high complexity, we propose a new and more efficient method, named the lattice approximation (LA) method.

A. Downlink SNR Distribution
For any given UAV location u and by dropping (u) for brevity, the SNR γ dl in (19) can be rewritten as which is mainly determined by the UAV-GBS association and the aggregate co-channel interference I. However, it is evident that they are coupled with each other. In the following, we first derive a general formula to obtain the downlink SNR cdf, by resolving the above coupling based on conditional probabilities.
From (22), the cdf of γ dl can be defined as where y > 0 is assumed; the right-hand side (RHS) of (a) takes expectation of the conditional probability P I ≥ −α 0 + C is /y|C is over the realization of the channel power gain C is with the associated GBS; and F I|Ci s (·) is the cdf of I conditioned on the realization of C is .
The realization of C is depends on the UAV-GBS association, whose pmf can be obtained similarly by Algorithm 1 in Section III. The overall algorithm to obtain the cdf of the downlink SNR γ dl is summarized in Algorithm 2. Denote T as the running time of each iteration (which will be specified later in Section IV-B) to obtain the conditional interference cdf F I|Ci s (·) in Step 4 of this algorithm. As a result, the worsecase complexity to compute the cdf of downlink SNR γ dl (u) at UAV location u is O(T |B|). Note that similar to Algorithm 1, we can stop the algorithm early if the current LoS channel power gain C L,im is smaller than C max NL , which is omitted in Algorithm 2 for brevity. The early stopping threshold for Algorithm 1 can also be similarly used to further reduce the number of iterations.  Obtain the conditional cdf F I|C is (·) by using the LA method in Section IV-B. 5: end for 6: Compute the cdf of γdl by (23).
The main challenge of implementing Algorithm 2 lies in how to obtain the conditional cdf F I|Ci s (·) in each iteration m, which is addressed in the following. For any given UAV location u and by dropping (u) for brevity, the aggregate interference power in (18) can be expressed as (24) which depends on two sets of independent RVs, i.e., the random channel activity µ i and the probabilistic LoS channel occurrence δ i of all the co-channel GBSs i ∈ B fs \ i s . Conditioned on the realization of channel power gain C is of the associated GBS i s , a subset of LoS channel realizations δ i , i ∈ B are implied. Specifically, for the m-th iteration in Algorithm 2, some of the co-channel GBSs should have the NLoS channel realization, i.e., δ ij = 0, ∀j = 1, · · · , m−1. As a result, the involved co-channel interference terms I ij , j = 1, · · · , m − 1, i j ∈ B fs \ i s , can only take two possible values between 0 (corresponding to µ i = 0) and C NL,ij (corresponding to µ i = 1) with probabilities 1 − ω dl,ij and ω dl,ij , respectively. On the other hand, the interference term I i from each remaining co-channel GBS i ∈ B fs \ i s can take three possible values of 0, C NL,i and C L,i with probabilities 1 − ω dl,i , ω dl,i (1 − p L,i ) and ω dl,i p L,i , respectively.
In general, consider a discrete RV ζ i which takes values from 1, · · · , L with probabilities p i,1 , · · · , p i,L , respectively, where L l=1 p i,l = 1. Based on independent ζ i , i = 1, · · · , M , a new RV is defined as which is named as the generalized Poisson multinomial (GPM) RV 10 , where z i L l=1 a i,l 1 l (ζ i ) can take L possible values a i,1 , · · · , a i,L , and 1 l (x) is the indicator function where 1 l (x) = 1 when x = l and 1 l (x) = 0 otherwise. Without loss of generality, assume that a i,1 ≤ a i,2 · · · ≤ a i,L . From (24) and (25), it follows that the conditional aggregate interference I can be modeled as a GPM variable.
For real-valued a i,l 's in practice, the value space of the composite GPM variable Z has a worst-case size of L M . As a result, in order to completely characterize the pmf of Z, an enumeration-based method has time complexity O(L M ), and moreover it requires additional O(L M ) time for sorting the mass points in order to obtain the cdf. To avoid such exponential complexity, in the sequel we propose the LA method to obtain the approximated cdf of Z efficiently with high accuracy.

B. Lattice Approximation (LA) Method
Instead of directly computing the pmf/cdf, we leverage the characteristic function (cf) 11 of Z, which always exists 10 It is called "generalized" as it allows each individual summand z i to take values from different real-valued sample spaces with different probabilities, thus extending the Poisson multinomial distribution [34] and further the ordinary multinomial distribution. 11 The cf of an RV X is defined as ϕ X (s) and uniquely characterizes its distribution. The cf approach is useful in the analysis of linear combinations of independent RVs. In our setting, the GPM variable Z is the sum of independent RVs z i , i = 1, · · · , M , whose cf is given by which is decomposed as the product of individual cf ϕ zi (s) given by After decomposition, computing the cf ϕ Z (s) at one point s has only linear complexity O(LM ), which involves computing M independent cfs ϕ zi (s), i = 1, · · · , M , using the given values of p i,l and a i,l for i = 1, · · · , M and l = 1, · · · , L. Next, we introduce a cf inversion lemma which will be used to recover the cdf of Z. Lemma 1. If an RV X has all its masses concentrated on the integer lattice N {0, 1, · · · , N − 1} of finite length N , then its pmf q n = P{X = n}, n ∈ N can be recovered by performing inverse discrete Fourier transform (DFT) on the N equally-spaced samples of the cf ϕ X (s), i.e., where Proof: If X has its masses concentrated on N with the pmf q n P{X = n}, n ∈ N , then its cf is given by which is continuous in s and has a period 2π (analogous to the frequency domain signal), and coincides with the discrete time Fourier transform (DTFT) for the "time" sequence q n , n ∈ N . To obtain perfect recovery (i.e., no aliasing) of the sequence q n , n ∈ N , we can take N samples from ϕ X (s) which are equally spaced around the unit circle as which is the DFT of the sequence q n , n ∈ N . Therefore, we can use the inverse DFT to recover the pmf sequence as in (28) and hence lemma 1 follows. Note that the fast Fourier transform (FFT) algorithm can be applied to perform the inverse DFT efficiently. In our context, the cf samples for the GPM variable Z can be obtained using (26). However, Lemma 1 cannot be applied directly to recover the pmf of Z. This is because the real-valued discrete GPM variable Z has its masses concentrated on a set of L M points, which are usually not equally spaced to form a lattice 12 or even integer lattice, and moreover the number of mass points to be recovered is overwhelming even for the FFT algorithm. To resolve the above challenge, we propose the LA method, where the original Z is approximated by an integer lattice with a bounded size of value space, and furthermore the efficient FFT algorithm is applied on its cf samples to recover its pmf and hence the approximated cdf of Z.
Specifically, we propose the offset-scale-quantize operation to convert a general non-lattice GPM variable Z into an integer-lattice GPM variableZ, and then use the cdf ofZ to obtain an approximated cdf of Z. Denote A 0 M i=1 a i,1 as the minimum possible value of Z, and A M i=1 a i,L −A 0 as the absolute range of possible values of Z. The first step is to offset each summand z i by a i,1 , respectively, so that each summand has a minimum value of 0 and so does their sum. This offsetting process is practically helpful for the stability of numerical computations involving complex numbers of potentially large phase (e.g., computing the cf samples in (26)), and also facilitates the subsequent scaling operation. Second, we perform a scale-and-quantize operation on z i − a i,1 and convert it intoz i , which is obtained by multiplying z i − a i,1 with a common factor of β > 0, and then rounding all the possible values β(a i,l − a i,1 )'s to their nearest integer values, denoted asã i,l , respectively, to which the original mass p i,l now attributes 13 . The absolute value range of the converted GPM variableZ M iz i can then be denoted as which is around βA with · denoting the ceiling function. As a result, the new GPM variableZ has its masses concentrated on the lattice N = {0, 1, · · · ,Ã} of finite length A + 1 ≈ βA + 1. Compared to the exponential number L M of mass points in a general non-lattice GPM variable Z, the scale-and-quantize operation effectively merges adjacent nonlattice mass points into an integer lattice with bounded size of value space, where the size can be controlled by the scaling factor β. Note that the offset-scale-quantize operation takes effect on each possible value of each summand z i , which has an overall complexity of O(LM ).
For the converted integer-valued GPM variableZ, samples of its cf ϕZ(s) can be obtained similar to (26). Therefore, we can apply Lemma 1 to obtain the exact pmf ofZ whose mass points are already sorted, and hence we can directly obtain its cdf FZ(·). The cdf of the original GPM variable Z can then be approximated as where (a) is due to possible quantization errors.
Note that a larger scaling factor β corresponds to more lattice points and hence smaller quantization errors. However, the time to compute βA samples of the cf ϕZ(s) in (26) increases in O(βALM ), while the time to perform FFT increases in O βA log 2 (βA) . Therefore, there is a trade-off between the accuracy and computational efficiency. The overall algorithm to obtain the cdf of a real-valued GPM variable is summarized in Algorithm 3, whose computational time T has the complexity of O(LM ) + O(βALM ) + O βA log 2 (βA) . Note that β can be chosen as β = c 0 /A so that βA is a large enough constant c 0 (e.g., in the range of 100 to 1000) to provide high approximation accuracy, and hence the LA method has a linear complexity in M .
Algorithm 3 LA method to approximate the cdf of non-lattice GPM Input: For each independent summand zi, i = 1, · · · , M , given its possible values a i,l with probabilities p i,l for l = 1, · · · , L. Let ai,1 ≤ ai,2 · · · ≤ ai,L. Output: The approximated cdf of the GPM variable Z = M i=1 zi.  [22] are used in the simulation for the channel model in (6).
A. Uplink Communication 1) Spatial Distribution of Non-Outage Probability: The non-outage probability 1 − p out,ul (u) at location u can be obtained by the analysis in Section III. The spatial distribution of uplink non-outage probability for the UAV located in the reference cell 0 at different altitude H u under different GBS antenna downtilt angle θ tilt (as in Fig. 3 (b) and (c)) is plotted in Fig. 4. It can be seen that the GBS antenna pattern, especially the sidelobes and nulls in between, has a significant impact on the spatial distribution of non-outage probability, where the UAV is mainly served by the sidelobe peaks of GBS 0 or other GBSs complementarily. Moreover, a larger |θ tilt | leads to overall smaller GBS antenna gain for the UAV above the GBS height, which results in overall lower strength of the communication link and hence lower coverage probability (average non-outage probability), as will be also shown later in Fig. 5. On the other hand, as the UAV altitude increases, the UAV-GBS link distance increases, while so does the LoS probability. This two factors affect the link strength and hence the coverage probability in opposite ways.
The spatial distribution of non-outage probability is helpful in facilitating UAV path planning/movement control. For the   Fig. 4(a), in order to move across the cell from location a to location c, the straight path a − c needs to fly over the outage region while the dashed red trajectory a−b−c enjoys full uplink coverage along the way with short traveling distance. Other examples may suggest UAV movement across different layers of altitude in the 3D space. Similar applications exist for the downlink case (e.g., based on Fig. 8 in the sequel).

example in
2) Coverage Probability: The uplink coverage probability p c,ul at an altitude H u can be obtained in (16). The trends of p c,ul versus the UAV altitude H u under different GBS antenna downtilt angle θ tilt and UAV antenna beamwidth Φ u are plotted in Fig. 5. For the cases with Φ u = 90 • (isotropic UAV antenna), the coverage probability p c,ul corresponds to the spatial average of non-outage probability in Fig. 4 and tends to decrease with the UAV altitude H u , while the oscillation is mainly due to the non-uniform GBS antenna pattern in the elevation domain. In contrast, for the cases with Φ u = 75 • , the UAV directional antenna gain is helpful for p c,ul (see Fig.  5 at high altitude, e.g., H u ≥ 80 m). On the other hand, at low altitude, the limited coverage range of the UAV antenna mainlobe limits its chance to be served by GBSs with good channel. As a result, the coverage probability first increases and then decreases with the UAV altitude.
In addition, the trend of p c,ul versus the SNR threshold η ul is plotted in Fig. 6. It can be seen that p c,ul in general decreases with η ul at a given UAV altitude. On the other hand, as the UAV altitude H u increases, the LoS probability also increases, resulting in an overall less variation of p c,ul versus η ul .
B. Downlink Communication 1) Aggregate Interference Distribution: Assume that the UAV is equipped with isotropic antenna of unit gain, located at horizontal location (0.3D, 0.1D) with altitude H u = 100 m, and is associated with the GBS with the largest LoS channel gain. We apply the LA method, enumeration method, Monte Carlo (MC) based simulation, and a benchmark Gaussian Approximation (GA) method to obtain the aggregate interference distribution under different cell loading factor ω dl , for comparison. For the MC method, 10 6 random samples of the aggregate interference I is generated in order to provide a good approximation for the true distribution, where each sample of I is drawn by summing over one realization of the interference terms I i , i ∈ B fs \ i s which are randomly and independently generated. On the other hand, the GA method is based on the central limit theorem, where the cdf of I is approximated by the non-negative part of the Gaussian cdf which is normalized such that the total probability is 1. The mean and standard deviation of I need to be computed in the GA method, based on the given values of a i,l 's and p i,l 's and with a computation time T of O(LM 2 ).
Consider a network area with radius D max = 3D, where there are 37 GBSs in B and 11 co-channel GBSs. The results are plotted in Fig. 7. It can be seen that the LA method matches almost exactly with the enumeration method, while both results are verified by MC simulation. On the other hand, the GA method provides fair approximation result for the case with moderate loading (e.g., ω dl = 0.5), while the approximation result is poor for the case with low (or high) loading. The average CPU time for the LA, GA, MC, and enumeration methods are 0.040, 0.028, 7.36 and 194 seconds, respectively, which are performed in MATLAB2015b on a laptop computer with Intel i7 2.7GHz CPU and 8GB memory without multi-core tasking. It can be seen that both the LA and GA methods run much faster than the MC and enumeration methods. Note that due to exponential time complexity, the enumeration method cannot be used for a larger setup with one or more tiers of co-channel GBSs (i.e., M ≥ 17). In summary, the LA method provides highly accurate approximation for the aggregate interference distribution with low time complexity, thus it will be applied in the rest of simulations.
2) Spatial Distribution of Non-Outage Probability: The downlink non-outage probability 1 − p out,dl (u) at location u can be obtained by the analysis in Section IV. The spatial distribution of downlink non-outage probability for the UAV located in the reference cell 0 at different altitude H u is plotted in Fig. 8. Compared to the uplink, the UAV's nonoutage probability in the downlink is determined by the direct communication link as well as the aggregate interference distribution. In particular, a larger downtilt angle |θ tilt | leads to overall smaller GBS antenna gain for the UAV above the GBS height, which reduces the overall strength of both direct link and interference links. However, the reduction on the aggregate interference is more significant, resulting in overall higher downlink coverage probability for the UAV above a certain altitude (e.g., H u ≥ 60 m), which is also shown next in Fig. 9.
3) Coverage Probability: The trends of downlink coverage probability p c,dl versus the UAV altitude H u under different UAV antenna beamwidth Φ u and GBS antenna downtilt angle θ tilt are plotted in Fig. 9. The curve oscillation is due to the non-uniform GBS antenna pattern which also results in complicated interference distribution. It can be seen that at high altitude, the directional UAV antenna is helpful for improving p c,dl , which not only brings directional antenna gain, but also effectively limits the interference from GBSs outside the UAV's antenna coverage range, hence resulting in higher p c,dl compared to the case with isotropic UAV antenna.
In addition, the trend of p c,dl versus the SNR threshold η dl is plotted in Fig. 10. It can be seen that at a given altitude, p c,dl decreases with η dl , and is lower for a higher loading factor ω dl .  Moreover, a higher loading factor ω dl also leads to overall less variation of p c,dl versus η dl .

VI. CONCLUSIONS
This paper studies the 3D system modeling and coverage performance analysis for network-connected UAVs in the cellular uplink and downlink communications. A 3D cellular network model is presented which incorporates the 3D airground channel and 3D patterns of the GBS/UAV antennas. Based on it, an analytical framework is further proposed to characterize the uplink/downlink 3D coverage performance, which effectively reduces the exponential complexity due to UAV-GBS association and coupled interference distribution. The conditional discrete interference is modeled as a new GPM RV, and a novel LA method is proposed to obtain the interference distribution efficiently with high accuracy. The 3D coverage analysis is validated by extensive numerical results, which also show the effects of cell loading factor, GBS antenna downtilt, UAV altitude and antenna beamwidth, and applications for UAV path planning/movement control. Our analytical framework is applicable to heterogeneous GBS locations, heights, antenna patterns and loading factors. The results based on directional UAV antenna and practically downtilted GBS antenna provide motivation for more advanced antenna and beamforming design, such as vertically sectorized GBS antenna and 3D digital beamforming at the GBS/UAV. Advanced interference mitigation techniques such as multi-cell coordinated GBS selection, channel allocation, power control and joint transmission/reception can also be applied to further improve the coverage performance, which will lead to promising future studies.