Downlink Coverage and Rate Analysis of an Aerial User in Vertical Heterogeneous Networks (VHetNets)

TL;DR

This paper analyzes the downlink coverage probability and rate for aerial users in vertical HetNets, considering aerial and terrestrial base stations, using stochastic geometry to derive exact and approximate expressions.

Contribution

It introduces a comprehensive stochastic geometry model for aerial users in VHetNets, incorporating LoS and NLoS channels, and derives simplified, accurate expressions for coverage and rate.

Findings

Aerial users at low altitudes connect to aerial-BSs in dense environments.

Directive beamforming improves aerial user performance by reducing interference.

Performance degrades in dense terrestrial networks despite beamforming.

Abstract

In this paper, we analyze the downlink coverage probability and rate of an aerial user in vertical HetNets (VHetNets) comprising aerial base stations (aerial-BSs) and terrestrial-BSs. The locations of terrestrial-BSs are modeled as an infinite 2-D Poisson Point Process (PPP), while the locations of aerial-BSs are modeled as a finite 2-D Binomial Point Process (BPP). Our cellular-to-air (C2A) channel model incorporates line-of-sight (LoS) and non-LoS transmissions between terrestrial-BSs and a typical aerial user, while we assume LoS transmissions for all aerial links. We assume that the aerial user is associated with an aerial-BS or terrestrial-BS that provides the strongest average received power. Using stochastic geometry, we derive exact and approximate expressions of the coverage probability and rate in terms of interference power's Laplace transform. The expressions are simplified assuming only LoS transmissions for the C2A channels. This enables easy-to-compute equations with good accuracy at elevated aerial user heights. We find that aerial users hovering at low altitudes tend to connect to aerial-BSs in denser terrestrial environments. Employing directive beamforming at aerial-BSs guarantees an acceptable performance at the aerial user by reducing interference signals received from the aerial-BSs. In denser terrestrial networks, the performance at the aerial user degrades substantially despite beamforming.

Tables

Table 1. TABLE I: The main relevant stochastic geometry-based works on VHetNets

Ref.	Aerial network	Terrestrial network	C2A links model	Coverage prob.	Rate	Aerial user
[14]	Finite 2-D BPP	Infinite 2-D PPP	LoS links	✓	✗	✓
[17]	Infinite 2-D PPP	N/A	LoS/NLoS links	✓	✓	✗
[18]	Infinite 2-D PPP	Infinite 2-D PPP	LoS/NLoS links	✓	✓	✗
[19]	Finite 2-D BPP	N/A	LoS links	✓	✗	✗
[23]	N/A	Infinite 2-D PPP	LoS/NLoS links	✓	✓	✓
[24]	N/A	Infinite 2-D PPP	LoS/NLoS links	✓	✗	✓
[29]	N/A	Infinite 2-D PPP	LoS/NLoS links	✓	✓	✓
[30]	Single aerial-BS	Single terrestrial-BS	LoS/NLoS links	✓	✗	✗
This paper	Finite 2-D PPP	Infinite 2-D PPP	LoS/NLoS links	✓	✓	✓

Table 2. TABLE II: Approximated LoS Probabilities for Different Environments

Environment	${𝐡}_T$ (meters)	a	b	c
Suburban	30	1	6.581	1
Urban	19	1	0.151	1
Dense urban	25	1	0.106	1
Highrise urban	62	1.124	0.049	1.024

Table 3. TABLE III: List of Symbols

Notation	Description
$h_{\mathrm{UT}}$	Height of the aerial user with reference to the terrestrial-BSs
$h_{\mathrm{UA}}$	Height of the aerial user with reference to the aerial-BSs
$h_U$	Height of the aerial user with reference to the ground level
$P_T,P_A$	Transmit power of terrestrial and aerial BSs, respectively
$r_D$	Radius of the circle where aerial-BSs are distributed
$H_C$	Height of the cylinder where aerial-BSs are distributed in the 3-D setup
${\Phi }_T,{\Phi }_A,{\Phi }_L,{\Phi }_N$	Tier of terrestrial-BSs, aerial-BSs, LoS terrestrial-BSs, or NLoS terrestrial-BSs, respectively
${\lambda }_T$	Density of terrestrial-BSs
$P_L(\cdot ),P_N(\cdot )$	Probability of the aerial user being in LoS and NLoS with terrestrial-BSs, respectively
$G_m^A,g_s^A$	Gain of the aerial-BS’s antenna mainlobe and sidelobe, respectively
$G_m^T,g_s^T$	Gain of the terrestrial-BS’s antenna mainlobe and sidelobe, respectively
$q_A$	Probability that the mainlobe of the interfering aerial-BS is directed toward the aerial user
$m_L,m_N,m_A$	Nakagami-$m$ fading parameter for LoS terrestrial-BS, NLoS terrestrial-BS, and aerial-BS, respectively
${\alpha }_L,{\alpha }_N,{\alpha }_A$	Path-loss exponent parameter for LoS terrestrial-BS, NLoS terrestrial-BS, and aerial-BS, respectively
$d_{L,x_i},d_{N,x_i},d_{A,x_i}$	Distance between the aerial user and an LoS terrestrial-BS, NLoS terrestrial-BS, or aerial-BS, respectively, located at $x_i$
$R_L,R_N,R_A$	Distance between the aerial user and its nearest LoS terrestrial-BS, NLoS terrestrial-BS, or aerial-BS, respectively
${\tilde{R}}_L,{\tilde{R}}_N,{\tilde{R}}_A$	Distance between the aerial user and its serving BS assuming the that the aerial user is associated with an LoS terrestrial-BS, NLoS terrestrial-BS, or aerial-BS, respectively
${ℰ}_L,{ℰ}_N,{ℰ}_A$	Event that describes an aerial user’s serving BS as LoS terrestrial-BS, NLoS terrestrial-BS, or aerial-BS, respectively
$H_L^{x_i},H_N^{x_i},H_A^{x_i}$	Channel power gain between the aerial user and an LoS terrestrial-BS, NLoS terrestrial-BS, or aerial-BS, respectively, located at $x_i$
$I_L,I_N,I_A$	Interference at the aerial user from LoS terrestrial-BS, NLoS terrestrial-BS, and aerial-BS, respectively
${𝒜}_L,{𝒜}_N,{𝒜}_A$	Probability of association between aerial user and LoS terrestrial-BS, NLoS terrestrial-BS, or aerial-BS, respectively

Table 4. TABLE IV: Interferer distances

Condition	${\tau }_{L|{ℰ}_{\nu }}\left(r\right)$	${\tau }_{N|{ℰ}_{\nu }}\left(r\right)$	${\tau }_{A|{ℰ}_{\nu }}\left(r\right)$	N’
${ℰ}_L$ ($x_0\in {\Phi }_L$)	$r$	where ${\zeta }_N^L={\left(\frac{{\eta }_L{h}_{\mathrm{UT}}^{{\alpha }_N}}{{\eta }_N}\right)}^{\frac{1}{{\alpha }_L}}$	where ${\zeta }_A^L\left(x\right)={\left(\frac{{\mu }_L{x}^{{\alpha }_A}}{{\mu }_A}\right)}^{\frac{1}{{\alpha }_L}}$	N
${ℰ}_N$ ($x_0\in {\Phi }_N$ )	${\left(\frac{{\mu }_L{r}^{{\alpha }_N}}{{\mu }_N}\right)}^{\frac{1}{{\alpha }_L}}$	$r$	$h_{\mathrm{UA}}$	N
${ℰ}_A$ ($x_0\in {\Phi }_A$)	where ${\zeta }_L^A={\left(\frac{{\mu }_A{h}_{\mathrm{UT}}^{{\alpha }_L}}{{\mu }_L}\right)}^{\frac{1}{{\alpha }_A}}$	$h_{\mathrm{UT}}$	$r$	N-1

Table 5. TABLE V: Parameter values

Parameter	Value	Parameter	Value	Parameter	value
$(P_T,P_A)$	(43, 30) dBm	($m_L,m_N,m_A$)	(2, 1, 2)	N	10
$r_D$	2000 meters	$(G_m^T,g_s^T)$	(0, -15) dB	(${\eta }_L,{\eta }_N,{\eta }_A$)	(-3, -20, -1) dB
$(G_m^A,g_s^A)$	(0, -20) dB	C2A	LoS/NLoS	${\sigma }^2$	-113 dBm
$({\alpha }_L,{\alpha }_N,{\alpha }_A)$	(2.5, 3.5, 2)	${\lambda }_T$	5	$T$	5 dB
Policy	N-OSS	Environment	Urban	${\theta }_B^A$	${18}^{\circ }$
($h_U,h_A$)	(50, 300) meters	$q_A$	0.1	$h_T$	20 meters

Equations

$$P_{\mathrm{L}}(z)=\prod _{n=0}^m\left[1-\exp \left(-\frac{{\left[h_{\mathrm{TX}}-\frac{\left(n+\frac{1}{2}\right)(h_{\mathrm{TX}}-h_{\mathrm{RX}})}{m+1}\right]}^2}{2{\delta }^2}\right)\right],$$

$$P_{\mathrm{L}}(\theta )=-a\exp \left(-b\theta \right)+c,$$

$$P_{\mathrm{L}}(z)=-a\exp \left(-b{\tan }^{-1}\left(\frac{h_{\mathrm{UT}}}{z}\right)\right)+c.$$

$$f_{{\mathrm{H}}_{\nu }^{x_j}}(x)=\frac{m_{\nu }^{m_{\nu }}{x}^{m_{\nu }x}}{\Gamma (m_{\nu })}{e}^{-m_{\nu }x},\nu \in \{\mathrm{L},\mathrm{N}\},$$

$$P_{r,j}^{\nu }=P_T{\eta }_{\nu }{G}^{\mathrm{T}}(r_j){\mathrm{H}}_{\nu }^{x_j}{d}_{\nu ,x_j}^{-{\alpha }_{\nu }},\nu \in \{\mathrm{L},\mathrm{N}\},$$

$$G^{\rm T}(r_{j})=\left\{\begin{array}[]{ll}G_{\rm m}^{\rm T},&r_{j}\in\mathcal{D}_{\rm TBS}\\ g_{\rm s}^{\rm T},&r_{j}\notin\mathcal{D}_{\rm TBS},\end{array}\right.$$

$$G^{\rm A}(d_{A,x_{i}})=\left\{\begin{array}[]{ll}G_{\rm m}^{\rm A},&-\frac{\theta_{\rm B}^{\rm A}}{2}\leqslant\psi\leqslant\frac{\theta_{\rm B}^{\rm A}}{2}\\ g_{\rm s}^{\rm A},&\text{otherwise},\end{array}\right.$$

$$P_{r,i}^A=P_A{G}^{\mathrm{A}}(d_{A,x_i}){\eta }_A{\mathrm{H}}_{\mathrm{A}}^{x_i}{d}_{A,x_i}^{-{\alpha }_A},$$

$$x_0=\text{arg}\max \{{\mu }_L{R}_L^{-{\alpha }_L},{\mu }_N{R}_N^{-{\alpha }_N},{\mu }_A{R}_A^{-{\alpha }_A}\},$$

$$\gamma=\left\{\begin{array}[]{ll}\frac{\mu_{L}{\rm H}_{\rm L}^{x_{0}}R_{L}^{-\alpha_{L}}}{I+\sigma^{2}},&{\rm if}\hskip 4.0ptx_{0}\in\Phi_{\rm L}(\mathcal{E}_{\rm L})\\ \frac{\mu_{N}{\rm H}_{\rm N}^{x_{0}}R_{N}^{-\alpha_{N}}}{I+\sigma^{2}},&{\rm if}\hskip 4.0ptx_{0}\in\Phi_{\rm N}(\mathcal{E}_{\rm N})\\ \frac{\mu_{A}{\rm H}_{\rm A}^{x_{0}}R_{A}^{-\alpha_{A}}}{I+\sigma^{2}},&{\rm if}\hskip 4.0ptx_{0}\in\Phi_{\rm A}(\mathcal{E}_{\rm A}),\end{array}\right.$$

$$I=\left\{\begin{array}[]{ll}I_{L}+I_{N}+I_{A},&\text{N-OSS}\\ I_{L}+I_{N},&\text{OSS},\hskip 8.0pt\text{ if }x_{0}\in\Phi_{\rm L}\cup\Phi_{\rm N}\\ I_{A},&\text{OSS},\hskip 8.0pt\text{ if }x_{0}\in\Phi_{\rm A},\end{array}\right.$$

$$I_L=\sum _{x_j\in {\Phi }_L\setminus x_0}{\mu }_L{\mathrm{H}}_{\mathrm{L}}^{x_j}{d}_{L,x_j}^{-{\alpha }_L},I_N=\sum _{x_j\in {\Phi }_N\setminus x_0}{\mu }_N{\mathrm{H}}_{\mathrm{N}}^{x_j}{d}_{N,x_j}^{-{\alpha }_N},\text{ and, }{I}_A=\sum _{i=1,x_i\in {\Phi }_A\setminus x_0}^N{P}_A{G}^{\mathrm{A}}(d_{A,x_i}){\eta }_A{\mathrm{H}}_{\mathrm{A}}^{x_i}{d}_{A,x_i}^{-{\alpha }_A},$$

$$f_{R_{\nu }}(r)$$

$$F_{R_{\nu }}(r)$$

$$f_{R_A}(r)=N\left(\frac{2r}{r_{\mathrm{D}}^2}\right){\left(\frac{d^2-r^2}{r_{\mathrm{D}}^2}\right)}^{N-1},h_{\mathrm{UA}}⩽r⩽d$$

$$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!F_{R_{A}}(r)=\left\{\begin{array}[]{ll}0,&r\leqslant h_{\rm UA}\\ 1-\left(\frac{d^{2}-r^{2}}{r_{\rm D}^{2}}\right)^{N},&h_{\rm UA}\leqslant r\leqslant d\\ 1,&r\geqslant d,\end{array}\right.$$

$${\mathcal{A}}_{\mathrm{L}}={\Xi }_{\mathrm{N}}^{\mathrm{L}}\times {\Xi }_{\mathrm{A}}^{\mathrm{L}},$$

$$\Xi^{\rm L}_{\rm N}=\left\{\begin{array}[]{ll}F_{R_{L}}(\zeta_{\rm N}^{\rm L})\!+\!\displaystyle\int_{\zeta_{\rm N}^{\rm L}}^{\infty}F_{R_{N}}^{(c)}\left(\!\!\left(\frac{\eta_{N}r^{\alpha_{L}}}{\eta_{L}}\right)^{\frac{1}{\alpha_{N}}}\!\!\right)f_{R_{L}}(r)dr,&h_{\rm UT}\leqslant\zeta_{\rm N}^{\rm L}\\ \displaystyle\int_{\zeta_{\rm N}^{\rm L}}^{\infty}F_{R_{N}}^{(c)}\left(\left(\frac{\eta_{N}r^{\alpha_{L}}}{\eta_{L}}\right)^{\frac{1}{\alpha_{N}}}\right)f_{R_{L}}(r)dr,&h_{\rm UT}\geqslant\zeta_{\rm N}^{\rm L},\end{array}\right.$$

$$\Xi^{\rm L}_{\rm A}=\left\{\begin{array}[]{ll}1+\displaystyle\int_{\zeta}^{\zeta_{\rm A}^{\rm L}(d)}\!\!\!\!F^{(c)}_{R_{A}}\!\left(\!\left(\!\frac{\mu_{A}r^{\alpha_{L}}}{\mu_{L}}\!\right)^{\frac{\!1}{\!\alpha_{A}}}\!\right)\!f_{R_{L}}\!(r)dr-F_{R_{L}}^{(c)}(\zeta),&h_{\rm UT}\!\leqslant\!\zeta_{\rm A}^{\rm L}(d)\\ 0,&h_{\rm UT}\geqslant\zeta_{\rm A}^{\rm L}(d),\end{array}\right.$$

$${\mathcal{A}}_{\mathrm{A}}={\Xi }_{\mathrm{N}}^{\mathrm{A}}\times {\Xi }_{\mathrm{L}}^{\mathrm{A}},$$

$$\Xi^{\rm A}_{\rm N}=\left\{\begin{array}[]{ll}\displaystyle\int_{h_{\rm UA}}^{d}\!F_{R_{N}}^{(c)}\left(\!\left(\!\frac{\mu_{N}r^{\alpha_{A}}}{\mu_{A}}\!\right)^{\frac{\!1}{\!\alpha_{N}}}\!\right)\!f_{R_{A}}\!(r)dr,&\zeta_{\rm N}^{\rm A}\leqslant\!h_{\rm UA}\\ F_{R_{A}}(\zeta_{\rm N}^{\rm A})\!+\!\!\displaystyle\int_{\zeta_{\rm N}^{\rm A}}^{d}\!\!\!F_{R_{N}}^{(c)}\left(\!\left(\!\frac{\mu_{N}r^{\alpha_{A}}}{\mu_{A}}\!\right)^{\frac{\!1}{\!\alpha_{N}}}\!\right)\!f_{R_{A}}\!(r)dr,&h_{\rm UA}\!\leqslant\!\zeta_{\rm N}^{\rm A}\!\leqslant\!d\\ 1,&\zeta_{\rm N}^{\rm A}\geqslant\!d,\end{array}\right.$$

$$\Xi^{\rm A}_{\rm L}\!=\!\left\{\begin{array}[]{ll}\displaystyle\int_{h_{\rm UA}}^{d}\!F_{R_{L}}^{(c)}\left(\!\left(\!\frac{\mu_{L}r^{\alpha_{A}}}{\mu_{A}}\!\right)^{\frac{\!1}{\!\alpha_{L}}}\!\right)\!f_{R_{A}}\!(r)dr,&\zeta_{\rm L}^{\rm A}\leqslant\!h_{\rm UA}\\ F_{R_{A}}(\zeta_{\rm L}^{\rm A})\!+\!\!\displaystyle\int_{\zeta_{\rm L}^{\rm A}}^{d}\!\!F_{R_{L}}^{(c)}\left(\!\left(\!\frac{\mu_{L}r^{\alpha_{A}}}{\mu_{A}}\!\right)^{\frac{\!1}{\!\alpha_{L}}}\!\right)\!f_{R_{A}}\!(r)dr,&h_{\rm UA}\leqslant\!\zeta_{\rm L}^{\rm A}\!\leqslant\!\!d\\ 1,&\zeta_{\rm L}^{\rm A}\geqslant\!d.\\ \end{array}\right.$$

$$\!\!\!f_{\tilde{R}_{L}|\mathcal{E}_{\rm L}}(r)\!\!=\!\!\left\{\!\!\!\begin{array}[]{ll}\frac{1}{\mathcal{A}_{\rm L}}F^{(c)}_{R_{N}}\left(\tau_{N|\mathcal{E}_{\rm L}}(r)\right)\times F^{(c)}_{R_{A}}\left(\tau_{A|\mathcal{E}_{\rm L}}(r)\right)\!\!\times\!\!f_{R_{L}}(\!r\!),&r\in\mathfrak{R}_{\rm L}\\ 0,&\text{otherwise},\end{array}\right.$$

$$\!\!\!f_{\tilde{R}_{A}|\mathcal{E}_{\rm A}}(r)\!\!=\!\!\left\{\!\!\!\begin{array}[]{ll}\frac{1}{\mathcal{A}_{\rm A}}F^{(c)}_{R_{N}}\left(h_{\rm UT}\right)\times F^{(c)}_{R_{L}}\left(\tau_{L|\mathcal{E}_{\rm A}}(r)\right)\times f_{R_{A}}(r),&r\in\mathfrak{R}_{\rm A}\\ 0,&\text{otherwise},\end{array}\right.$$

$$f_{\tilde{R}_{N}|\mathcal{E}_{\rm N}}(r)=\left\{\!\!\!\begin{array}[]{ll}\frac{1}{\mathcal{A}_{\rm N}}F^{(c)}_{R_{L}}\left(\tau_{L|\mathcal{E}_{\rm N}}(r)\right)\times F^{(c)}_{R_{A}}\left(h_{\rm UA}\right)\times f_{R_{N}}(r),&r\in\mathfrak{R}_{\rm N}\\ 0,&\text{otherwise},\end{array}\right.$$

$$\mathcal{L}_{(I_{N}+I_{L})|\mathcal{E}_{\nu}}(s)\!=\!\!\!\!\!\!\prod_{\omega\in\{\rm L,N\}}\!\!\!\!\!\exp\Big{(}\!\!-\!2\pi\lambda_{\rm T}\int_{\tau_{\omega|\mathcal{E}_{\rm\nu}}(r)}^{\infty}\!\!\Big{(}\!1\!-\!\left(\!\frac{\!\!\!m_{\omega}}{m_{\omega}\!+\!s\mu_{\omega}t^{-\alpha_{\omega}}}\right)^{\!\!\!m_{\omega}}\!\Big{)}tP_{\rm\omega}\!\!\left(\!\!\sqrt{\!t^{2}\!-\!h_{\rm UT}^{2}}\right)\!dt\!\Big{)}.$$

$${\mathcal{L}}_{I_A|{\mathcal{E}}_{\nu }}(s)$$

$$\Omega(x,g)=\frac{x^{\alpha_{A}m_{A}+2}}{g^{m_{A}}}{\rm G}_{2,2}^{1,2}\Biggl{[}\!\frac{m_{A}x^{\alpha_{A}}}{sP_{A}g\eta_{A}}\!\Bigg{|}\!\ \ {\!1\!-\!m_{A}\!-\!\frac{2}{\alpha_{A}},1\!-\!m_{A}\!\atop 0,-m_{A}-\frac{2}{\alpha_{A}}}\Biggr{]},$$

$$f_{R_L}(r)$$

$$F_{R_L}(r)$$

$$\mathcal{A}_{\rm L}\!=\!\!\left\{\begin{array}[]{ll}1\!\!-\exp\left(-\pi\lambda_{\rm T}\left(\zeta^{2}-h_{\rm UT}^{2}\right)\right)+\frac{\pi\lambda_{\rm T}h_{\rm UT}^{2}}{r_{\rm D}^{2N}}\displaystyle\sum_{i=0}^{N}\binom{N}{i}(-1)^{i}d^{2(N-i)}\left(\frac{\mu_{A}}{\mu_{L}}\right)^{\frac{2i}{\alpha_{A}}}\left(\pi\lambda_{\rm T}\right)^{-\frac{\alpha_{L}i}{\alpha_{A}}}\\ \times\left[\Upsilon(\zeta)-\Upsilon(\zeta_{\rm A}^{\rm L}(d))\right],&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!h_{\rm UT}\!\leqslant\!\zeta_{\rm A}^{\rm L}(d)\\ 0,&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!h_{\rm UT}\geqslant\zeta_{\rm A}^{\rm L}(d).\end{array}\right.$$

Full text

I Introduction

Unmanned aerial vehicles (UAVs, also referred to as drones) have drawn a great deal of interest in recent years due to their flexibility and cost effectiveness in executing a range of applications, including search-and-rescue missions, aerial imaging, surveillance, and package-delivery, to name few [1]. To accommodate such high-capacity and demanding scenarios, ubiquitous coverage for UAVs is becoming crucial for delivering ultra-reliable and low-latency real-time communications [2, 3]. On this point, it may be worth mentioning that 3GPP in its recent report [4] discussed the great potential of UAVs for providing cellular connectivity in 5G networks. More specifically, the report investigated the communication needs of UAVs (referred to as on-board radio access nodes (UxNB) in [4]). Connected with this, a new paradigm has recently emerged where existing cellular networks (e.g., Long Term Evolution (LTE) networks) and future cellular networks (5G and beyond) can be re-utilized to provide cellular connectivity to UAVs (aerial users) [5]. Indeed, some field trials have been carried out to substantiate the feasibility of terrestrial networks for supporting aerial user communications [6]. The authors of [7], for instance, provided a comprehensive study of the potential of terrestrial base stations (terrestrial-BSs) and aerial base stations (aerial-BSs) in providing cellular coverage to aerial users based on 4G deployment. Besides the multitude of applications where UAVs are deployed as aerial users, UAVs can also be equipped with mounted antennas (i.e., as aerial-BSs) to provide cellular connectivity in case of disasters where the terrestrial network is in failure [8, 9]. Also, significantly, aerial-BSs can be deployed to support terrestrial network connectivity for high-capacity demands, such as concerts and sporting events. Due to the nature of temporary spikes in throughput demand, deploying aerial-BSs is the fastest and most cost-effective solution, compared to terrestrial-BSs [10, 11, 12, 13, 14]. Aerial-BSs have also been investigated as a promising solution for providing pervasive cellular connectivity in remote areas where the installation of conventional terrestrial network infrastructure is deemed challenging [15]. Clearly, the deployment of aerial users is imminent due to their flexibility and the great improvements of their payload capacity and extended flight duration. However, the performance of aerial users in terms of coverage and achievable throughput have not received much attention in the literature on integrated aerial and terrestrial networks (i.e., vertical heterogeneous networks or VHetNets) [14, 16, 17, 18]. In this paper, we leverage tools from stochastic geometry to derive closed-form expressions for the downlink coverage and achievable rate of a typical aerial user in VHetNets under the assumptions of a 2-D infinite terrestrial network (terrestrial-BSs distributed over an infinite 2-D region with a certain density) and a 2-D finite aerial network. The results provide some interesting insights for the future deployment of aerial users/BSs.

I-A Motivation and Related Work

Due to the emerging applications of UAVs as aerial-BSs and aerial users, the performance of VHetNets involving aerial and terrestrial users have attracted a great deal of interest. While the research community has extensively studied the efficiency of deploying UAVs as aerial-BSs [19, 20, 18, 21], the study of aerial users in terms of coverage probability and achievable rates has been lacking [6, 22, 23, 24]. For instance, the authors in [6] presented a non-orthogonal multiple access (NOMA) technique for uplink communications between aerial users and terrestrial-BSs. In [25], the authors investigated interference-aware path planning for aerial users using a deep reinforcement learning algorithm. Moreover, the work in [6] tackled the interference problem in aerial communications. A thorough comparison between the existing network infrastructure and futuristic deployment scenarios involving massive multiple-input multiple-output (MIMO) systems on the performance of aerial users was studied in [22]. While both [6] and [22] focused on terrestrial-BSs as servers for aerial users, neither explored the opportunity of deploying aerial-BSs for the cellular connectivity of aerial users.

Based on the tools of stochastic geometry, the authors in [23] and [24] investigated the feasibility of serving aerial and terrestrial users with a terrestrial wireless network. The authors claimed that the line-of-sight (LoS) condition between aerial users and terrestrial-BSs degrades the signal-to-interference-and-noise ratio SINR at the terrestrial users due to the strong interference signals received from the LoS aerial users. The same authors also presented several tuning parameters to combat these strong interference signals. However, they did not consider the possibility of using aerial-BSs capable of serving aerial users in the event of a congested terrestrial wireless network.

While the Poisson point process (PPP) assumption has become the baseline for modeling the locations of terrestrial-BSs, it is not suitable for modeling aerial networks, especially those with a small number of aerial-BSs deployed over a limited area. Since the Binomial point process (BPP) was adopted to model ad hoc terrestrial networks with a given number of nodes distributed over a circular area [26, 27], it presents a reasonable model for aerial networks while also assuring the tractability of mathematical analysis. It is noteworthy that the deployment of BSs (terrestrial or aerial) is not random in practice. However, by assuming BPP aerial network, the mathematical analysis remains tractable. For instance, the work in [19] was the first to consider a BPP model for aerial networks. In this work, the links between terrestrial users and aerial-BSs were assumed to be in LoS condition (non-LoS (NLoS) links were not considered). Such an assumption may not be applicable to aerial-BSs operating at very low heights and/or in dense environments due to the high likelihood of NLoS occurrences in these scenarios. In addition, the work in [19] limited the association of terrestrial users to only aerial-BSs, and it did not consider the possibility of the terrestrial user being served by a nearby terrestrial-BS that would provide a better channel condition than that provided by an aerial-BS, if there were one. Similar to the authors in [24], the authors in [28] investigated the opportunity of reusing existing terrestrial networks to serve aerial users. Their work suggested several interference mitigation techniques both in the uplink and downlink transmissions in order to maintain an acceptable performance of the terrestrial network for both terrestrial and aerial users. In [29], the use of coordinated multi-point transmissions to provide seamless cellular connectivity for aerial users was studied. In the same work, clustered small cell BSs were considered for serving aerial users in content-caching architecture, where popularly requested content was cached for aerial users which could then be transmitted to terrestrial users. However, the work was limited to scenarios where terrestrial-BSs did not coexist with aerial-BSs. The authors in [17] proposed a framework to analyze the coverage and rate of a VHetNet. Yet the analysis was exclusively for a typical terrestrial user; so it is not applicable for aerial users.

In Table I, we summarize the stochastic geometry-based works in VHetNets in terms of characteristics for terrestrial and aerial networks, C2A transmissions link type, derived performance, and whether they account for aerial users. It may be worth mentioning that aerial-BSs can serve terrestrial users along with their aerial counterparts. However, in this paper, we focus on filling the gap in the literature by investigating only the aerial user’s performance in a VHetNet setup.

I-B Contributions

The contributions of this paper can be summarized as follows:

• •

We propose a novel framework to analyze the downlink coverage and rate of an aerial user served by a VHetNet comprising of aerial-BSs and terrestrial-BSs. The model consists of terrestrial-BSs that follow an infinite 2-D PPP and aerial-BSs that follow a finite 2-D BPP deployed at a particular height above the ground [19]. We model the cellular-to-air (C2A) links according to the general framework provided by the International Telecommunications Union (ITU) in its recommendation report [31], which incorporates LoS and NLoS transmissions occurring according to a given probability. By contrast, the air-to-air (A2A) channels are assumed to be in LoS condition due to the absence of obstacles between the aerial user and aerial-BSs. In order to alleviate the strong interference signals received from the aerial-BSs, we incorporate directional beamforming at the aerial-BSs, where the interfering main beam of the aerial-BS is directed toward the typical aerial user with a certain probability depending on its antenna beamwidth. Two spectrum sharing policies between aerial-BSs and terrestrial-BSs (referred to as orthogonal spectrum sharing (OSS) and non-OSS (N-OSS)) are considered.

• •

Since the expression of the LoS probability of C2A links is difficult to use, and the widely used air-to-ground (A2G) model in [20] and [32] was exclusively developed for terrestrial users, we derive an easy-to-use compact exponential function as an approximation using curve fitting. We subsequently show that the obtained simple LoS probability matches perfectly with the ITU model in [31].

• •

We derive exact and approximate expressions for the coverage probability and rate in terms of the Laplace transform of interference power, and we show that the derived expressions match the simulations perfectly.

• •

Under the assumption that there are only LoS C2A links (no NLoS transmissions), closed-form expressions for the aerial user association probabilities and the Laplace transform of the terrestrial interference power are derived.

• •

The mathematical analysis of the interference received from the aerial-BSs is challenging, especially because of the beamforming that directs the main beam of the aerial-BS toward the typical aerial user with a certain probability. This assumption results in two dependent tiers of aerial-BSs, where the number of aerial-BSs in each tier depends on the other tier. Based on this setup, we derive an exact expression for the Laplace transform of the interference received from aerial-BSs in terms of the Meijer-G function. To the best of our knowledge, there are no existing studies that present a closed-form expression for the interference received from aerial-BSs. In fact, closed-form expressions are often presented as integral expressions [19]. In this case, the derived closed-form expression simplifies the evaluation of the expressions of the coverage probability and rate.

This paper is organized as follows. The system model is presented in Section II. In Section III, we derive expressions for the association probabilities and the Laplace transforms of the aggregated interference powers. Also in Section III, we derive some useful distance distributions that are used in subsequent sections. Closed-form expressions for the coverage probability and the average achievable rate are derived in Section IV. Finally, we validate and compare the derived analytical expressions with Monte Carlo simulations and investigate the impact of several system parameters on the performance of the VHetNets.

II System Model

We assume a network of single tier terrestrial-BSs, denoted by $\Phi_{\rm T}$, uniformly distributed at ground level (i.e., terrestrial-BSs follow PPP) with density $\lambda_{\rm T}$ [terrestrial-BSs/km2] and height $h_{\rm T}$ 111Although terrestrial-BSs are often a mix of macro-, micro-, and pico-BSs, we only consider a single tier terrestrial network composed of macro-BS since macro-BSs transmit at higher powers than other terrestrial-BSs and are thus more suitable for providing wireless connectivity to aerial users.. Since in practice, aerial-BSs can hover at different heights, we start by assuming that a finite number of aerial-BSs, $N$, hovers at different heights confined between minimum and maximum heights, $h_{\rm min}$ and $h_{\rm max}$, respectively, as shown in Fig. 1(a). We also assume that $N$ aerial-BSs are uniformly distributed within a cylinder with a radius $r_{\rm D}$ and a height $\rm H_{C}$ to form a BPP with a center $(0,0,\frac{h_{\rm min}+h_{\rm max}}{2})$. Introducing such a randomness to the height of aerial-BSs makes the distance distributions between a typical aerial user and any aerial-BS difficult to derive without much of a useful outcome. Since it has been shown in [19] that the performance of a 3-D BPP aerial network matches perfectly that of 2-D BPP aerial-BSs deployed at an average height, $h_{\rm A}=\frac{h_{\rm min}+h_{\rm max}}{2}$, we assume in our framework, $N$ aerial-BSs hovering at a height $h_{\rm A}$ from the ground and uniformly distributed over a disc with a center $(0,0,h_{\rm A})$ and radius $r_{\rm D}$, as illustrated in Fig. 1(b). So although we place the aerial-BSs at the same height $h_{\rm A}$, our framework is still applicable for aerial-BSs with different heights. Simulation results will substantiate the correctness of this assumption. It may be worth mentioning that the deployment of BSs (either terrestrial of aerial) is not random in practice. However, the point process assumption enables the tractability of the analysis of the network using stochastic geometry tools while tracking realistic BS deployment in practice [33, 17, 19]. In this work, our focus is mainly on an aerial user’s performance in VHetNets. The transmit power of the terrestrial-BSs and aerial-BSs are assumed to be $P_{T}$ and $P_{A}$, respectively. Finally, we consider a typical aerial user hovering at a height $h_{\rm U}$, where $h_{\rm T}<h_{\rm U}<h_{\rm A}$.

II-A Cellular-to-air channel

The widely used A2G channel model proposed in [32] and [20] is not applicable to aerial users because it was developed for terrestrial users positioned at heights of 1.5 meters, while our system model involves terrestrial-BSs positioned at much greater heights (e.g., 30 meters for macro-BSs) [34]. Therefore, in this section, we develop a more tractable expression for the LoS probability. According to the ITU recommendation report [31], the probability of an LoS between a transmitter and receiver with heights $h_{\rm TX}$ and $h_{\rm RX}$, respectively, is given by [31]

[TABLE]

where $z$ denotes the horizontal distance between the transmitter and the receiver on the $xy$ plane, and $m=\lfloor\left(\frac{z\sqrt{\alpha\beta}}{1000}-1\right)\rfloor$ where $\alpha$, $\beta$, and $\delta$ are environment-related coefficients given in Table I in [35].

As we can see in (1), the LoS probability $P_{\rm L}(z)$ is not a continuous function of $z$ which makes the analysis intractable. Similar to the work in [32] and [20], we simplify (1) to a more tractable exponential function under some assumptions. This can be done by approximating the curve of $P_{\rm L}(z)$ in (1) by a continuous exponential function of the elevation angle $\theta$ for a fixed terrestrial-BS height. From Fig. 1(b), the horizontal distance $z$ can be written as $z=(h_{\rm U}-h_{\rm T})/\tan(\theta)$. Comparing the parameters in (1) with the parameters of the system model presented in Fig. 1(b), we have, $h_{\rm TX}=h_{\rm T}$ and $h_{\rm RX}=h_{\rm U}$. Notably, the resulting plot of the expression in (1) will smooth for large values of $h_{\rm U}$ (we used $h_{\rm U}=10,000$ meters in our simulation). Thus, the probability of an LoS condition becomes a continuous function of the elevation angle $\theta$ and the environment parameters. In Fig. 2, we plotted the ITU LoS probability for four different environments with a solid red line. As we can see, the trend of the curves can be closely approximated to a continuous exponential function. Using the Matlab® feature ‘Curve Fitting’, we have

[TABLE]

where $\theta=\tan^{-1}\left(\frac{h_{\rm U}-h_{\rm T}}{z}\right)$, and $a$, $b$, and $c$ are parameters that depend on the environment and height of the terrestrial-BSs. The LoS probability model used in this work is presented in Table II. Indeed, the LoS probability derived in (2) for C2A links presents a compact expression that depends only on the environment and the height of the aerial user in its elevation angle for a given terrestrial-BS height.

Fig. 2 shows the LoS probability using the original ITU model in (1) and the approximated expression in (2) for different environments. As we can see, the approximated LoS expression in (2) matches the original ITU model very closely. Each C2A link between the aerial user, which is at a height $h_{\rm UT}=h_{\rm U}-h_{\rm T}$ from the terrestrial-BS, is assumed to be either an LoS or NLoS link with an LoS probability $P_{\rm L}(z)$ given by

[TABLE]

Finally, the NLoS probability is given by $P_{\rm N}(z)=1-P_{\rm L}(z)$.

Due to the fact that the aerial user is either in an LoS or NLoS condition with each terrestrial-BS [18], the set of terrestrial-BSs can be broken down into two independent inhomogeneous PPPs, where the LoS terrestrial-BSs form a subset $\Phi_{\rm L}$ with density $\lambda_{\rm T}P_{\rm L}(z)$, while the NLoS terrestrial-BSs form a subset $\Phi_{\rm N}$ with density $\lambda_{\rm T}P_{\rm N}(z)$ [36]. We assume that the LoS and NLoS C2A channels experience Nakagami-$m$ fading with different $m$ parameters, and therefore the received power from the terrestrial-BS located at point $x_{j}$ is Gamma-distributed, i.e., ${\rm H}_{\nu}^{x_{j}}\sim\text{Gamma}\left(m_{\nu},\frac{1}{m_{\nu}}\right)$, where $m_{\nu}$, $\nu\in\{\rm L,N\}$ denotes the fading parameters for the LoS and NLoS C2A links. The probability density function (PDF) of the Gamma-distributed channel gain is given by [18]

[TABLE]

where $\Gamma(\cdot)$ is the Gamma function defined as $\Gamma(m_{\nu})=\int_{0}^{\infty}t^{m_{\nu}-1}e^{-t}dt$.

The received power at the typical aerial user from a terrestrial-BS located at $x_{j}$ is given by

[TABLE]

where $\eta_{\nu}$ are the excess losses for the LoS and NLoS C2A links, respectively [32]. Moreover, $d_{\nu,x_{j}}$ denotes the distance between the typical aerial user and a terrestrial-BS from the tier $\Phi_{\nu},\nu\in\{\rm L,N\}$ located at $x_{j}$. $\alpha_{\nu}$, $\nu\in\{\rm L,N\}$, is the path-loss coefficient. Finally, $G^{\rm T}(r_{j})$ denotes the antenna directivity gain of the terrestrial-BS’s antenna located at $x_{j}$ and the typical aerial user, which can be written as [24]

[TABLE]

where $\mathcal{D}_{\rm TBS}$ denotes the set of all distances $r_{j}$ satisfying $r_{j}\tan\left(\theta_{t}+\frac{\theta_{\rm B}^{\rm T}}{2}\right)<h_{\rm U}<h_{\rm T}-r_{j}\tan\left(\theta_{t}-\frac{\theta_{\rm B}^{\rm T}}{2}\right)$, where $r_{j}$ refers to the horizontal distance separating the aerial user from the projection of the terrestrial-BS location on the $xy$ plane, and $\theta_{t}$ and $\theta_{\rm B}^{\rm T}$ are the down-tilt and beamwidth angles of the terrestrial-BS’s antenna. In other words, $G^{\rm T}(r_{j})$ determines whether the aerial user falls within the mainlobe or sidelobe of the terrestrial-BS’s antenna.

Since the terrestrial-BS antennas are typically tilted toward the ground in order to provide cellular coverage to terrestrial users [37], aerial users that hover at a greater height than the terrestrial-BS, i.e., $h_{\rm T}<h_{\rm U}$, will most likely receive signals from the terrestrial-BS antenna sidelobes. For the tractability of the analysis in this paper, we assume that $G^{\rm T}(r_{j})=g_{\rm s}^{\rm T},\ \forall\ r_{j}$. In the simulation results section, we will provide findings to substantiate the accuracy of this assumption.

II-B Air-to-air channel

Since the aerial user is assumed to hover above rooftops, the A2A links between the typical aerial user and the aerial-BSs are assumed to be in LoS condition [38, 39]. We also assume that the aerial-BSs employ directive beamforming to improve the SINR. As a result, the gain of the aerial-BS antenna located at $x_{i}$, denoted by $G^{\rm A}$, at the typical aerial user is given by [40]

[TABLE]

where $d_{A,x_{i}}$ denotes the distance between the typical aerial user and the aerial-BS located at $x_{i}$, $G_{\rm m}^{\rm A}$ and $g_{\rm s}^{\rm A}$ are the gains of the mainlobe and sidelobe, respectively, $\psi$ is the sector angle, and $\theta_{\rm B}^{\rm A}\in[0,180]$ is the beamwidth in degrees [41, Fig. 3]. Thus, the received power at the aerial user from an aerial-BS located at $x_{i}$ is given by

[TABLE]

where $\eta_{A}$ represents the excess losses, ${\rm H}_{\rm A}^{x_{i}}$ is the Gamma-distributed channel power gain, i.e., ${\rm H}_{\rm A}^{x_{i}}\sim\text{Gamma}\left(m_{A},\frac{1}{m_{A}}\right)$, with a fading parameter $m_{A}$, and $\alpha_{A}$ is the path-loss exponent. Unlike terrestrial-BSs whose antennas with fixed antenna patterns are deployed to serve terrestrial users, the aerial-BSs with dwon-tilted antennas can direct their beams toward their associated aerial users. Hence, the beamforming gain from the serving aerial-BS is always $G_{\rm m}^{\rm A}$ [41]. For the remaining interfering aerial-BSs, their main beams are not necessarily aligned with the typical aerial user since their beam orientations are uniformly distributed in $[0,180^{\circ}]$. We introduce a probability $q_{\rm A}$ depending on the beamwidth $\theta_{\rm B}^{\rm A}$ that quantifies the likelihood of an interfering aerial-BS’s mainlobe being directed toward the typical aerial user. It follows that $q_{\rm A}=\mathbb{P}\left(G^{\rm A}(d_{A,x_{i}})=G_{\rm m}^{\rm A}\right)=\frac{\theta_{\rm B}^{\rm A}}{180}$. The interfering aerial-BS gain is $g^{\rm A}_{\rm s}$ when the typical receiver is in its sidelobe direction with a probability of occurrence of $1-q_{\rm A}$ [41, Fig. 3]. Typically, $g^{\rm A}_{\rm s}$ is $20$ dB less than the mainlobe gain $G_{\rm m}^{\rm A}$ [41]. It is worth noting that realistic 3-D antenna radiation models, such as those presented in [42] and [43], can be adopted in our VHetNet’s framework analysis. However, as shown above, we choose to investigate a more analytically tractable BS antenna model. Future works can focus on more realistic 3-D antenna patterns in the context of VHetNets for even more practical insights.

The typical aerial user is assumed to be served by the BS (located at point $x_{0}$) that provides the strongest long-term averaged received power [44]. Indeed, since the A2A links are in LoS, it is possible for a distantly located aerial-BS to offer better SINR than that of a closer terrestrial-BS due to differences in path-loss parameters. It should be noted that if the aerial user is associated with a specific tier of BS, i.e., $\{\Phi_{\rm L},\Phi_{\rm N},\Phi_{\rm A}\}$, the serving BS would be the nearest BS to that specific tier. Thus, by assuming that the average power of all channels is 1, i.e., $\mathbb{E}[{\rm H}_{\rm L}^{x_{j}}]=\mathbb{E}[{\rm H}_{\rm N}^{x_{j}}]=\mathbb{E}[{\rm H}_{\rm A}^{x_{j}}]=1$ for each $x_{j}\in\Phi_{\rm T}\cup\Phi_{\rm A}$, the serving BS is given by

[TABLE]

where $\mu_{\nu}=P_{T}\eta_{\nu}g_{\rm s}^{\rm T}$, $\nu\in\{\rm L,N\}$, $\mu_{A}=P_{A}G_{\rm m}^{\rm A}\eta_{A}$, $R_{L}=\underset{\forall x_{j}\in\Phi_{\rm L}}{\min}d_{L,x_{j}}$, and, $R_{N}=\underset{\forall x_{j}\in\Phi_{\rm L}}{\min}d_{N,x_{j}}$, and $R_{A}=\underset{\forall x_{i}\in\Phi_{\rm A}}{\min}d_{A,x_{i}}$. As a result, the SINR at the typical aerial user is given by

[TABLE]

where $\mathcal{E}_{\nu}$, $\nu\in\{\rm L,N,A\}$, denotes the event that the serving BS belongs to the tier $\Phi_{\nu}$ and $\sigma^{2}$ is the additive white Gaussian noise power. Finally, $I$ refers to the aggregate interference power.

Two spectrum sharing policies between aerial-BSs and terrestrial-BSs are used: orthogonal spectrum sharing (OSS) and non-orthogonal spectrum sharing (N-OSS) [30]. OSS/N-OSS implies that aerial-BSs and terrestrial-BSs operate on the same/different frequencies. The interference $I$ is given by

[TABLE]

where

[TABLE]

where $G^{\rm A}(d_{A,x_{i}})=G_{\rm m}^{\rm A}$ with a probability of $q_{\rm A}$ and $G^{\rm A}(d_{A,x_{i}})=g_{\rm s}^{\rm A}$ with a probability of $1-q_{\rm A}$.

We provide the list of symbols in Table III.

III Relevant Distance Distributions and Association Probabilities

In this section, we begin by deriving several relevant distance distributions that will be useful in deriving the association probabilities and Laplace transforms of the interference powers given in (11). Afterwards, simplified expressions for the association probabilities and the terrestrial interference’s Laplace transforms are presented under the assumption that all C2A links are in LoS condition.

III-A Distance distributions of the nearest BSs

In Lemma 1, we present the distribution of the distances between the aerial user and its nearest aerial-BS, LoS terrestrial-BS, and NLoS terrestrial-BS.

Lemma 1**.**

The PDF and the cumulative distribution function (CDF) of $R_{\nu},\nu\in\{\rm L,N\}$ are given by

[TABLE]

Proof: See Appendix A.

The PDF and the CDF of $R_{A}$ are given by [19, eq. (7)]

[TABLE]

[TABLE]

where $d=\sqrt{r_{\rm D}^{2}+h_{\rm UA}^{2}}$.

III-B Distances of the nearest interfering BSs

In Table IV, we summarize the distances $\tau_{\nu|\mathcal{E}_{\nu}}(r),\nu\in\{\rm L,N,A\}$ between the typical aerial user and the nearest interfering BSs from the three tiers $\Phi_{\rm L},\Phi_{\rm N},\Phi_{\rm A}$ where the serving BS is at a distance of $r$. More specifically, the set of distances where $r\geqslant\tau_{\nu|\mathcal{E}_{\nu}}(r)$ quantify the region over which the interfering BSs from $\Phi_{\nu}$ exist. Note that $N^{\prime}$ in Table IV quantifies the number of interfering aerial-BSs.

III-C Association probabilities

As stated previously, the typical aerial user is associated with an aerial-BS, LoS terrestrial-BS, or NLoS terrestrial-BS with probabilities given in the following lemmas.

Lemma 2**.**

The probability that the typical aerial user is associated with an LoS terrestrial-BS is given by

[TABLE]

where

[TABLE]

and,

[TABLE]

with \zeta=\text{arg}\max\Bigl{\{}h_{\rm UT},\zeta_{\rm A}^{\rm L}(h_{\rm UA})\Bigr{\}}.

Proof: See Appendix B.

Lemma 3**.**

The probability of the typical aerial user being associated with an aerial-BS is given by

[TABLE]

where

[TABLE]

where $\zeta_{\rm N}^{\rm A}=\left(\frac{\mu_{A}h_{\rm UT}^{\alpha_{N}}}{\mu_{N}}\right)^{\frac{1}{\alpha_{A}}}$, and,

[TABLE]

Proof: The proof follows the same steps as that of Lemma 2; therefore, it is omitted here.

Finally, the probability of the aerial user being associated with an NLoS terrestrial-BS is given by $\mathcal{A}_{\rm N}=1-\mathcal{A}_{\rm L}-\mathcal{A}_{\rm A}$.

III-D Conditional distance distributions of the serving BS

In this section, we present the distribution of the distances between the typical aerial user and its serving BS, where the latter is an LoS terrestrial-BS, NLoS terrestrial-BS, or aerial-BS.

Lemma 4**.**

Where the aerial user is associated with an LoS terrestrial-BS (i.e., the event $\mathcal{E}_{\rm L}$ occurs), the distribution of the distance $\tilde{R}_{L}$ between the aerial user and the serving LoS terrestrial-BS is given by

[TABLE]

where $\mathfrak{R}_{\rm L}=[h_{\rm UT},\zeta_{\rm A}^{\rm L}(d)]$. Note that according to (19), $\mathcal{A}_{\rm L}=0$ outside $\mathfrak{R}_{\rm L}$.

Proof: See Appendix C.

Lemma 5**.**

Where the aerial user is associated with an aerial-BS (i.e., the event $\mathcal{E}_{\rm A}$ occurs), the distribution of the distance $\tilde{R}_{A}$ between the aerial user and the serving aerial-BS is given by

[TABLE]

where $\mathfrak{R}_{\rm A}=[h_{\rm UA},d]$.

Proof: The proof follows the same steps as that of Lemma 4; therefore, it is omitted here.

Lemma 6**.**

Where the aerial user is associated with an NLoS terrestrial-BS (i.e., the event $\mathcal{E}_{\rm N}$ occurs), the distribution of the distance $\tilde{R}_{N}$ between the aerial user and the serving NLoS terrestrial-BS is given by

[TABLE]

where $\mathfrak{R}_{\rm N}=[h_{\rm UT},d]$.

III-E Laplace transform of the aggregated interference

Depending on the spectrum sharing policy, the aerial user may receive interference signals from aerial-BSs and terrestrial-BSs whose Laplace transforms are characterized in the following lemmas.

Lemma 7**.**

Conditioned on the event $\mathcal{E}_{\nu}$, $\nu\in\{\rm L,N,A\}$, the Laplace transform of the terrestrial interference power, $I_{L}+I_{N}$, is given by

[TABLE]

Proof: See Appendix D.

Lemma 8**.**

Conditioned on the event $\mathcal{E}_{\nu}$, $\nu\in\{\rm L,N,A\}$, the Laplace transform of the aerial interference power, $I_{A}$, is given by

[TABLE]

where $\Omega[\cdot,\cdot]$ is given below by

[TABLE]

where ${\rm G}[\cdot]$ is the Meijer-G function given in [45, eq (9.301)].

Proof: See Appendix E.

III-F Simplified C2A channel model

The previous results (association probabilities and Laplace transform of interference) require numerical evaluations of multiple integrals. These expressions can be simplified by noting that the aerial user may hover at elevated heights resulting in LoS transmissions with the terrestrial-BSs. Therefore, we assume that the C2A links are in LoS condition with the aerial user. The validity of this assumption will be investigated with simulations in Section V. Under this assumption, the PDF and CDF of the distance between the aerial user and the nearest LoS terrestrial-BS are given by [41, eqs. (8) and (9)]

[TABLE]

Corollary 1**.**

With the assumption that the aerial user is in LoS with the terrestrial-BSs, the association probability with a terrestrial-BS given in (17) can be written as in (31) at the top of the next page, where $\Upsilon(x)=\Gamma\left[\frac{\alpha_{L}}{\alpha_{A}}i+1,\pi\lambda_{\rm T}x^{2}\right]$, with $\Gamma[\cdot,\cdot]$ denoting the upper incomplete Gamma function [45, eq (8.358.2)]. Thus, the probability of association with an aerial-BS is given by $\mathcal{A}_{\rm A}=1-\mathcal{A}_{\rm L}$.

Proof: Substitute (29) into (19) and apply the Binomial theorem. Then, using the definition of the upper incomplete Gamma function and some manipulations, we obtain (31).

Corollary 2**.**

The Laplace transform of the terrestrial interference, $\mathcal{L}_{I_{L}|\mathcal{E}_{\nu}}(s)$, where $\nu\in\{\rm L,A\}$, is given by

[TABLE]

where $\Theta(\cdot,\cdot)$ is defined as

[TABLE]

Proof: Starting from (26), substituting $P_{L}(t)=1,\quad t\geqslant h_{\rm UT}$ and applying the identity $a^{n}-b^{n}=(a-b)\sum_{i=0}^{n-1}a^{i}b^{n-1-i}$, [46, eq (1.43)], and [45, eq (7.811.2)] yield (32).

IV Coverage and Rate Analysis

IV-A Downlink coverage probability

The coverage probability is defined as the probability of the SINR at the typical aerial user exceeding a predetermined threshold $T$. By applying the law of total probability, the coverage probability at the typical aerial user is given by

[TABLE]

where $\mathcal{A}_{\nu}$ and $\mathcal{C}_{\nu},\nu\in\{\rm L,N,A\}$ denote the association probability with the $\Phi_{\nu}$ tier and the coverage probability given that the aerial user is associated with the BS from $\Phi_{\nu}$, respectively.

Theorem 1**.**

The coverage probability $\mathcal{C}_{\nu}$ conditioned on the events $\mathcal{E}_{\nu},\nu\in\{\rm L,N,A\}$ is given by

[TABLE]

where $V=I+\sigma^{2}$ and $\mathcal{L}_{V|\mathcal{E}_{\nu}}(s)$ is given by $\mathcal{L}_{V|\mathcal{E}_{\nu}}(s)=\exp(-\sigma^{2}s)\mathcal{L}_{I|\mathcal{E}_{\nu}}(s)$.

Note that under the N-OSS policy, $\mathcal{L}_{I|\mathcal{E}_{\nu}}(s)=\mathcal{L}_{(I_{N}+I_{L})|\mathcal{E}_{\nu}}(s)\times\mathcal{L}_{I_{A|\mathcal{E}_{\nu}}}(s)$.

Proof: See Appendix F.

IV-B Downlink average achievable rate

In this section, we derive the average achievable rate of the typical aerial user. Following the same analysis as the coverage probability, the average ergodic rate is given by

[TABLE]

where $\mathcal{R}_{\nu}$ is the average achievable rate of the aerial user when it is associated with a BS from the tier $\Phi_{\nu}$.

Theorem 2**.**

The average achievable rate of the typical aerial user when it is associated with $\Phi_{\nu}$ tier is given by

[TABLE]

Proof: Using Shannon’s theorem, the average achievable rate can be expressed as

[TABLE]

where (a) follows from $\mathbb{E}[Z]=\int_{0}^{\infty}\mathbb{P}(Z>t)dt$. Substituting (10) into (38) and following the same steps as those of Theorem 1 yields (37).

IV-C Downlink coverage and average achievable rate approximations

The evaluation of the results presented in (35) and (37) for the conditional coverage probabilities and average achievable rates involves higher order derivative of the Laplace transform, which leads to a more complex computation especially for large fading parameters, i.e., $m_{\nu}$, $\nu\in\{\rm L,N,A\}$. Thus, a lower bound of the Gamma distribution’s CCDF is used instead of the exact evaluation.

Theorem 3**.**

The approximation of the conditional coverage probability $\mathcal{C}_{\nu}^{\rm approx}$, $\nu\in\{\rm L,N,A\}$ is given by

[TABLE]

where $\rho_{\nu}=m_{\nu}(m_{\nu}!)^{\frac{-1}{m_{\nu}}}$, with $m_{\nu}>1$.

Proof: Using [47], the CCDF of a Gamma-distributed random variable, $\rm H_{\nu}^{x_{0}}$, can be tightly lower bounded by

[TABLE]

where the equality holds when $m_{\nu}=1$. Thus, substituting this upper bound to (c) in (49) along the same mathematical manipulations concludes the proof.

Theorem 4**.**

The conditional average achievable rate approximation $\mathcal{R}_{\nu}^{\rm approx}$, $\nu\in\{\rm L,N,A\}$ is given by

[TABLE]

Proof: The proof follows the same steps as that of Theorem 3; therefore, it is omitted here.

V Results and Discussion

In this section, we use Monte Carlo simulations to validate the analytical expressions. We also investigate the impact of several design parameters on the VHetNet’s performance. Table V summarizes the parameters used in the simulations, unless otherwise indicated. It should be mentioned that we used parameters set by the 3GPP report in [48] which limit the maximum hovering height for aerial users to 300 meters. The height of an aerial-BS, however, can extend up to a few kilometers [19, 20, 17].

In Fig. 3, we compare aerial-BS and terrestrial-BS association probabilities with different aerial user heights and environments. Notably, at any aerial user height, $\mathcal{A}_{\rm A}+(\mathcal{A}_{\rm L}+\mathcal{A}_{\rm N})=1$. At low aerial user heights, the aerial user tends to be associated with an aerial-BS because the C2A links are dominated by NLoS links, which have poor channel conditions in comparison to LoS A2A links. As the height increases, the terrestrial-BS association probability increases since more terrestrial-BSs come within LoS condition of the aerial user. However, with a further increase in height, the terrestrial-BS association probability starts to decrease because of the significant path-loss caused by the increasing distances between the aerial user and the terrestrial-BSs. Essentially, for a given environment and height, the aerial user probability of association to either a terrestrial or aerial network indicates the likelihood of the serving BS belonging to either of the two networks, respectively.

Fig. 4 shows the coverage probability versus the SINR threshold for different environments. For a particular SINR threshold, as the environment becomes less dense, more terrestrial-BSs come into LoS with the aerial user (i.e., more interfering LoS terrestrial-BSs) which increases the interference. As discussed in Section II, the analytical tractability of a VHetNet with aerial-BSs at varying heights is quite challenging. Interestingly, aerial-BSs hovering at the same height $h_{\rm A}$ guarantee a comparable performance to VHetNets involving aerial-BSs with heights uniformly distributed between $h_{\rm max}$ and $h_{\rm min}$ with an average height of $h_{\rm A}$. The results presented in Fig. 4 show that a small variation in the height of aerial-BSs does not affect the overall performance. This is because the aerial-BSs are distributed over the area of a disc with a radius that is much larger than the aerial-BS heights, which means that small variations in aerial-BS heights do not have a significant impact on the distribution of the distances separating the aerial user from the aerial-BSs.

Fig. 5 shows the performance of VHetNets where aerial-BSs hover at different heights but are confined between two values: $h_{\rm min}=h_{\rm A}-\frac{\rm H_{C}}{2}$ and $h_{\rm max}=h_{\rm A}+\frac{\rm H_{C}}{2}$ as depicted in Fig. 1(a), where aerial-BSs are uniformly distributed within a finite cylinder of height $\rm H_{C}$ and radius $r_{\rm D}$. Fig. 5 shows that the coverage probability using a 3-D BPP aerial network is the same as its counterpart assuming a 2-D BPP at a fixed height $h_{\rm A}$ (in this case $\rm H_{C}=0$) when $\rm H_{C}$ is much larger than the aerial-BSs height, $h_{\rm A}$. However, the 2-D BPP assumption does not guarantee a likewise coverage probability under 3-D BPP when the radius over which the aerial-BSs are distributed uniformly, $r_{\rm D}$, is comparable to their heights, $h_{\rm A}$. Indeed, works in [19] and [49] indicate that aerial-BSs are distributed over a large area with a radius in the order of a few kilometers, whereas aerial-BS heights are in the order of a few hundred meters. Given this, we wish to emphasize that the framework presented in this work is applicable under the assumption that the radius of the disk $r_{\rm D}$ over which the aerial-BSs can be distributed is much larger than the cylinder height $\rm H_{\rm C}$. That is, an aerial network with a dense deployment of aerial-BSs, where their heights and spacing are comparable, does not necessarily exhibit the same performance as the aerial network in this work. This is because the assumption concerning order of magnitude differences between $r_{\rm D}$ and $\rm H_{C}$ is violated, and, therefore, revisiting the performance analysis using more comprehensive models is required. Indeed, under the assumption $r_{\rm D}\gg\rm H_{C}$, Fig. 5 shows that the complexity from assuming a 3-D BPP network can be avoided by adopting a 2-D BPP network with a fixed height without compromising the performance accuracy or correctness. Finally, we can affirm that the assumption of fixed heights of aerial-BSs is legitimate enough to enable the tractability of the analysis.

Fig. 6 illustrates the coverage probability for different terrestrial-BS tilt-angles. The performance of VHetNets that assume different $\theta_{t}$ is close enough to the performance of the system, assuming that the aerial user is always serviced by the terrestrial-BS antenna sidelobes $G^{\rm T}(r_{j})=g_{\rm s}^{\rm T}$, $\forall\ r_{j}$. This observation illustrates well the assumption that terrestrial-BSs cover aerial users mainly with their antenna sidelobes, since these are usually down-tilted toward terrestrial users. Further, as we pointed out, this assumption is the key enabler of the tractability of the association probabilities analysis.

In Fig. 7, we compare the coverage probability with different aerial user heights with OSS and N-OSS policies. Under the N-OSS policy, as the aerial user height increases, more terrestrial-BSs come within LoS of the aerial user, which increases the terrestrial interference at a greater rate than the increase in the desired signal power. Under the OSS policy, a different trend can be observed at low aerial user heights, i.e., the coverage probability increases. This is because at very low aerial user heights, the improvement in the desired signal power becomes greater than the increase in the terrestrial interference (since the aerial user may be in LoS with few terrestrial-BSs while terrestrial interference is dominated by NLoS terrestrial-BSs). However, with a further increase in the aerial user height, the interference received from the LoS terrestrial-BSs dominates, which degrades the coverage probability. Indeed, the significant difference in coverage probability under N-OSS and OSS policies is due to the heavy interference that the aerial user experiences under the N-OSS policy which includes aerial and terrestrial interferers alike. Under the OSS policy, by contrast, the aerial user is subject to interference only from BSs that belong to the same network (terrestrial or aerial) as its serving BS. It should also be mentioned that under the assumption that the C2A links are all LoS links, the coverage probability converges to the actual coverage probability.

Fig. 8 shows the coverage probability versus the density of the terrestrial-BSs for different aerial user heights and beamwidths under the OSS and N-OSS policies. Under the N-OSS policy, we notice that the coverage probability decreases as the density of terrestrial-BSs increases due to the increase in the terrestrial interference (more LoS terrestrial-BSs exist as the aerial user height increases). We can also observe that for a given density of terrestrial-BSs and regardless of the $q_{\rm A}$, the coverage probability decreases as the aerial user height increases due to the increase in the terrestrial interference (terrestrial interference is dominated by LoS terrestrial-BSs). By contrast, under the OSS policy and at a high aerial user height (250 meters), the density of the terrestrial-BSs has no impact on the coverage probability, since the aerial user is always associated with an aerial-BS according to Fig. 3, and thus, it does only receive interference from aerial-BSs. However, a similar trend is not observed if the aerial user is close to the terrestrial-BSs. For instance, at a low aerial user height (50 meters), the coverage probability decreases as the density of the terrestrial-BSs increases, since the aerial user may be associated with a terrestrial-BS, and so increasing the density of terrestrial-BSs increases the terrestrial interference. Overall, aerial BSs with more directed beams (smaller beamwidth implies smaller $q_{\rm A}$) improves the coverage probability because it decreases the aerial interference power.

We plot the coverage probability versus the number of aerial-BSs for different aerial user heights in Fig. 9. The approximation of the coverage probability using the tight lower bound of the Gamma distribution’s CCDF is in a good agreement with the exact performance. Generally, increasing the number of aerial-BSs decreases the coverage probability due to the fact that the increase in the aerial interference is greater than the increase in the desired signal. However, at high aerial user heights (e.g., 250 meters) and low aerial-BS numbers, the coverage probability increases as the number of aerial BSs increases. This is because the increase in aerial interference power is low (small number of aerial-BSs) compared to the improvement in the desired aerial signal. Overall, denser aerial-BSs decreases the coverage probability especially at low aerial user heights.

The validation of the LoS C2A assumption, which is presented in Section III-F, is shown in Fig. 10. At low aerial user heights (e.g., $h_{\rm U}=30$ meters), LoS C2A assumption is not a good approximation. This is because of the high likelihood of obstructions (NLoS links) between the aerial user and terrestrial-BSs. However, at elevated heights (e.g., $h_{\rm U}=270$ meters), the aerial user is very likely to be in LoS with terrestrial-BSs (no or a few NLoS terrestrial-BSs exist) which justifies the correctness and accuracy of the LoS C2A assumption.

In Fig. 11, we show the accuracy of the average achievable rate approximation given in (41). NLoS links are assumed to follow Rayleigh fading, $m_{N}=1$, whereas LoS links follow Nakagami-$m$ with $m_{L}=2$. We illustrate the average achievable rate for different fading parameter $m_{A}$ for the A2A channels. According to the figure, the approximate rate derived in (41) matches very closely the exact analytical results for a smaller fading parameter $m_{\rm A}$. As the latter increases, the approximation in (41) becomes an upper bound for the exact average achievable rate. Moreover, as the aerial-BSs height increase, the rate decreases. This can be explained by the fact that the power signal received from the aerial-BSs is weakened by the high path-loss, whereas the aggregate interference power received from the terrestrial-BSs remains the same.

VI Conclusion

Using stochastic geometry tools, we proposed a complete framework to analyze the coverage probability and rate of a typical aerial user served by a network of aerial-BSs and LoS/NLoS terrestrial-BSs. Useful distance distributions between the aerial user and its nearest terrestrial-BSs (LoS and NLoS) and aerial-BS were derived, which enabled the computation of the association probabilities of the typical aerial to terrestrial and aerial networks. Exact and approximations of the coverage probability and average achievable rate were presented assuming Nakagami-$m$ fading for all C2A and A2A channels. Simulations revealed that the approximations of the coverage probability and average achievable rate agree very well with the exact expressions, especially with a small fading parameter. Simple expressions were presented under the assumption of only LoS transmissions in the C2A links. Several conclusions can be drawn from the simulations and evaluated analytical expressions:

• •

The association probabilities of an aerial user with an aerial-BS or terrestrial-BS depend on the height of the aerial user, the height of the aerial-BS, and the environment. For instance, the aerial user tends to connect to an aerial-BS when the former is hovering at a very low height and in a dense environment. It was also shown that an optimum aerial user height exists that maximizes the association probability to a terrestrial-BS. As the aerial user’s height increases, the terrestrial-BS is no longer able to provide enough signal power due to high path-loss, which means the aerial user will then certainly be served by an aerial-BS.

• •

Simulation results show that the aerial user receives strong interference signals from LoS terrestrial-BSs and aerial-BSs, which degrade the coverage probability and rate. Orthogonal spectrum sharing between terrestrial and aerial networks considerably increases the coverage and rate performance of the aerial user. The simulations shows that an optimum height exists for which the coverage probability of the aerial user is maximized.

• •

It was also shown that employing directive beamforming at LoS aerial-BSs is essential in the operation of VHetNets as it decreases the aggregate power of aerial-BSs interference, and as a result it increases the aerial user coverage substantially. Other interference mitigation techniques might also be deployed to improve the aerial user operation.

• •

The results further show that aerial users at elevated heights are mainly in LoS with the terrestrial-BSs (no or a few NLoS links exist). Thus, the analysis can be simplified by assuming two tiers of BSs (i.e., LoS terrestrial-BSs and aerial-BSs). The simulation proves that this approximation is accurate at elevated heights and matches well the performance with LoS/NLoS C2A links. However, this approximation is not valid for low aerial user heights where NLoS links are more likely to exist and cannot be dismissed in the performance evaluation.

• •

Finally, it can be concluded from the results that increasing the density of the aerial-BSs may or may not degrade the performance, depending on the height of both the aerial-BSs and the aerial user. As we saw, a denser aerial network decreases the coverage probability of the aerial user, especially when it hovers at low heights.

Appendix A Proof of Lemma 1

$R_{\nu}$ can be written as $R_{\nu}=\sqrt{Z_{\nu}^{2}-h_{\rm UT}^{2}}$, $\nu\in\{\rm L,N\}$, where $Z_{\nu}$ is the horizontal distance between the aerial user and the nearest LoS and NLoS terrestrial-BSs. Using the null probability of PPP [33], we have

[TABLE]

Since $r=\sqrt{z^{2}+h_{\rm UT}^{2}}$, taking the integral in (a) with respect to $z$ completes the proof of (14) (note that $t$ is a dummy variable). (13) can be obtained by differentiating $F_{R_{\nu}}(r)$ with respect to $r$.

Appendix B Proof of Lemma 2

The association probability of an aerial user with an LoS terrestrial-BS can be written as

[TABLE]

where (a) follows from the independence of the two point processes that represent the aerial-BSs and terrestrial-BSs, (b) follows from the definition of the Complementary CDF and averaging over $R_{L}$. Finally, using (14) and (16) along with some mathematical manipulations completes the proof.

Appendix C Proof of Lemma 3

The CDF of ${\tilde{R}_{L}}$ can be written as

[TABLE]

where (a) follows from Bayes’ rule, (b) is obtained from (43), and (c) follows from averaging over $R_{L}$. Finally, differentiating (45) with respect to $r$ completes the proof.

Appendix D Proof of Lemma 7

The Laplace transform of the aggregated interference power is given by

[TABLE]

where (a) follows from (11) and the independence of the small-scale fading and PPP, and (b) is obtained from the moment generating function of the Gamma distribution. Finally, using the probability generating functional (PGFL) and the results in Table IV complete the proof.

Appendix E Proof of Lemma 8

It was proven in [19, Corollary 3] that when the serving BS (either terrestrial-BS or aerial-BS) is located at a distance $r$ from the aerial user, the distribution of the distance between the aerial user and the j-th interfering aerial-BS $d_{A,x_{j}},j\in[1,N^{\prime}]$ is given by

[TABLE]

where $N^{\prime}$ is the number of interfering aerial-BSs. It should be noted that there are $N^{\prime}-i$ interfering aerial-BSs that have a gain of $G_{\rm m}^{\rm A}$ with a probability $q_{\rm A}$ and $i$ interfering aerial-BSs that have a gain of $g_{\rm s}^{\rm A}$ with a probability of $1-q_{\rm A}$. Hence, $i$ follows a Binomial distribution $\mathcal{B}(N^{\prime},q_{\rm A})$. Indeed, the two subsets of interfering aerial-BSs (those with a gain of $G_{\rm m}^{\rm A}$ and those with a gain $g_{\rm s}^{\rm A}$) are dependent where the joint distribution of $N^{\prime}-i$ and $i$ is a multinomial distribution on $N^{\prime}$ trials with success probabilities $q_{\rm A}$ and $1-q_{\rm A}$, respectively.

The Laplace transform of the interference received from aerial-BSs is given by

[TABLE]

where (a) follows from using the Gamma distribution of the small scale fading channel gain of the A2A link. (b) follows from the binomial distribution of $i$ with $f_{i}(i)=\sum_{i=0}^{N^{\prime}}\binom{N^{\prime}}{i}q_{\rm A}^{N^{\prime}-i}(1-q_{\rm A})^{i}$. (c) is obtained after averaging over $d_{j}$. Using the identity $(1+z)^{a}=1/\Gamma(-a){\rm G}_{1,1}^{1,1}[z|{1-a\atop 0}]$ and [45, eq(7.811.2)] along with some additional mathematical manipulations, (27) is obtained.

Appendix F Proof of Theorem 1

The coverage probability conditioned on the event that the aerial user is connected to a BS from $\Phi_{\nu}$, $\nu\in\{\rm L,N,A\}$, is given by

[TABLE]

where (a) follows from averaging the coverage probability over the random variables, $\{{\rm H}_{\nu}^{x_{0}},\tilde{R}_{\nu},V\}$ and (b) is obtained by exploiting the independence between the three random variables. Moreover in (49), (c) follows after applying the CCDF of Gamma-distributed channel gain ${\rm H}_{\nu}^{x_{0}}$, and (d) is obtained by assuming $m_{\nu}$ to be an integer and using the series expansion of the upper incomplete Gamma function. Finally, using the identity, i.e., $\mathbb{E}_{V}\left[V^{k}\exp(-sV)\right]=(-1)^{k}\frac{\partial\mathcal{L}_{V}(s)}{\partial s^{k}}$, and averaging over $\tilde{R}_{\nu}$, we obtain (35).