Statistical Modeling of the FSO Fronthaul Channel for UAV-based Communications

TL;DR

This paper develops a statistical channel model for UAV-based free space optics (FSO) communication, accounting for UAV motion, orientation fluctuations, and non-orthogonal laser incidence, to improve system performance analysis and design.

Contribution

It introduces a novel channel model for UAV-based FSO systems that incorporates geometric and misalignment losses due to UAV motion and orientation fluctuations.

Findings

The model accurately predicts outage probability and ergodic rate.

Adjusting beam width can minimize outage probability under given fluctuations.

Simulations validate the analytical models and offer design insights.

Abstract

In this paper, we investigate the statistics of the free space optics (FSO) communication channel between a hovering unmanned aerial vehicle (UAV) and a central unit. Two unique characteristics make UAV-based FSO systems significantly different from conventional FSO systems with static transceivers. First, for UAV-based FSO systems, the incident laser beam is not always orthogonal to the receiver lens plane. Second, both position and orientation of the UAV fluctuate over time due to dynamic wind load, inherent random air fluctuations in the atmosphere around the UAV, and internal vibrations of the UAV. On the contrary, for conventional FSO systems, the laser beam is always perpendicular to the receiver lens plane and the relative movement of the transceivers is limited. In this paper, we develop a novel channel model for UAV-based FSO systems by quantifying the corresponding geometric and misalignment losses (GML), while taking into account the non-orthogonality of the laser beam and the random fluctuations of the position and orientation of the UAV. In particular, for diverse weather conditions, we propose different fluctuation models for the position and orientation of the UAV and derive corresponding statistical models for the GML. We further analyze the performance of a UAV-based FSO link in terms of outage probability and ergodic rate and simplify the resulting analytical expressions for the high signal-to-noise ratio (SNR) regime. Finally, simulations validate the accuracy of the presented analysis and provide important insights for system design. For instance, we show that for a given variance of the fluctuations, the beam width should be properly adjusted to minimize the outage probability.

Equations

$$y_s=h{x}_s+n,$$

$$h=\eta h_p{h}_a{h}_g,$$

$$h_p=10^{-\kappa L/10},$$

$$\alpha ={\left[\exp \left(\frac{0.49{\sigma }_R^2}{{\left(1+1.11{\sigma }_R^{12/5}\right)}^{7/6}}\right)-1\right]}^{-1},$$

$$\beta ={\left[\exp \left(\frac{0.51{\sigma }_R^2}{{\left(1+0.69{\sigma }_R^{12/5}\right)}^{5/6}}\right)-1\right]}^{-1}.$$

$$(x,y,z)=\mathbf{r}+ȷ\mathbf{d},$$

$$\displaystyle\mathbf{d}=\big{(}\sin\phi\cos\theta,\sin\phi\sin\theta,\cos\phi\big{)}.$$

$$\displaystyle\mathbf{b}=\Big{(}0,r_{y}-r_{x}\tan\theta,r_{z}-r_{x}\frac{\cot\phi}{\cos\theta}\Big{)}.$$

$$I^{\mathrm{orth}}(l;L)=\frac{2}{\pi w_L^2}\exp \left(-\frac{2l^2}{w_L^2}\right),$$

$$w_L=w_0\sqrt{1+\left(1+\frac{2w_0^2}{{\rho }^2(L)}\right){\left(\frac{\lambda L}{\pi w_0^2}\right)}^2}.$$

$$I(y,z)$$

$$h_g(\mathbf{r},\mathbf{\omega })=\underset{(y,z)\in \mathcal{A}}{\iint }I(y,z)\mathrm{d}y\mathrm{d}z,$$

$$h_g^{\mathrm{low}}(\mathbf{r},\mathbf{\omega })=\frac{2\sin \psi }{\pi w_L^2}\underset{(y,z)\in \mathcal{A}}{\iint }$$

$$\displaystyle\exp\bigg{(}-\frac{2}{w^{2}_{L}}\Big{(}(y-u)^{2}+\sin^{2}\phi\cos^{2}\theta{z}^{2}\Big{)}\bigg{)}\mathrm{d}y\mathrm{d}z$$

$$h_g^{\mathrm{upp}}(\mathbf{r},\mathbf{\omega })=\frac{2\sin \psi }{\pi w_L^2}\underset{(y,z)\in \mathcal{A}}{\iint }$$

$$\displaystyle\exp\bigg{(}-\frac{2}{w^{2}_{L}}\Big{(}\sin^{2}\phi\cos^{2}\theta(y-u)^{2}+{z}^{2}\Big{)}\bigg{)}\mathrm{d}y\mathrm{d}z,\,\,\quad$$

$${\tilde{h}}_g^{\mathrm{low}}(\mathbf{r},\mathbf{\omega })=A_0\exp \left(-\frac{2u^2}{t_1{w}_L^2}\right)$$

$${\tilde{h}}_g^{\mathrm{upp}}(\mathbf{r},\mathbf{\omega })=A_0\exp \left(-\frac{2u^2}{t_2{w}_L^2}\right),$$

$$A_0=\mathrm{erf}({\nu }_1)\mathrm{erf}({\nu }_2).$$

$$h_g(\mathbf{r},\mathbf{\omega })\approx A_0\exp \left(-\frac{2u^2}{t{w}_L^2}\right),$$

$$\mathbf{r}={\mathbf{\mu }}_{\mathbf{r}}+{\mathbf{ϵ}}_{\mathbf{r}}\text{and}\mathbf{\omega }={\mathbf{\mu }}_{\mathbf{\omega }}+{\mathbf{ϵ}}_{\mathbf{\omega }}.$$

$$\displaystyle\mu_{\theta}=\begin{cases}\pi+\tan^{-1}\big{(}\frac{\mu_{y}}{\mu_{x}}\big{)}&\mathrm{if}\,\,\mu_{x}>0\\ \tan^{-1}\big{(}\frac{\mu_{y}}{\mu_{x}}\big{)}&\mathrm{otherwise},\end{cases}\quad\text{and}$$

$$\displaystyle\mu_{\phi}=\pi-\cos^{-1}\bigg{(}\frac{\mu_{z}}{\sqrt{\mu_{x}^{2}+\mu_{y}^{2}+\mu_{z}^{2}}}\bigg{)}.$$

$${\mathbf{ϵ}}_{\mathbf{r}}={\mathbf{ϵ}}_{\mathbf{r}}^{\mathrm{IG}}\sim \mathcal{N}(0,{\mathbf{Q}}_{\mathbf{r}}^{\mathrm{IG}})\text{and}{\mathbf{ϵ}}_{\mathbf{\omega }}={\mathbf{ϵ}}_{\mathbf{\omega }}^{\mathrm{IG}}\sim \mathcal{N}(0,{\mathbf{Q}}_{\mathbf{\omega }}^{\mathrm{IG}}),$$

$$f_{h_g}(h)=\frac{\sqrt{t{w}_L^2}}{2\sqrt{2}h\sqrt{\ln \left(\frac{A_0}{h}\right)}}{f}_u\left(\sqrt{\frac{t{w}_L^2}{2}\ln \left(\frac{A_0}{h}\right)}\right).$$

$$f_{h_g}(h)=$$

$${\mathbf{\Sigma }}_{\mathrm{IG}}=\left[\begin{array}{cc}{\sigma }_y^2+c_1^2{\sigma }_x^2+c_2^2{\sigma }_{\theta }^2& c_1{c}_5{\sigma }_x^2+c_2{c}_4{\sigma }_{\theta }^2\\ c_1{c}_5{\sigma }_x^2+c_2{c}_4{\sigma }_{\theta }^2& {\sigma }_z^2+c_3^2{\sigma }_{\varphi }^2+c_4^2{\sigma }_{\theta }^2+c_5^2{\sigma }_x^2\end{array}\right].$$

$$\underset{h\to 0}{\lim }{f}_{h_g}(h)\stackrel{(a)}{=}$$

$${\mathbf{\Sigma }}_{\mathrm{IG}}=\left[\begin{array}{cc}{\sigma }_y^2+{\mu }_x^2{\sigma }_{\theta }^2& 0\\ 0& {\sigma }_z^2+{\mu }_x^2{\sigma }_{\varphi }^2\end{array}\right],$$

$$f_{h_g}(h)=$$

$${\mathbf{ϵ}}_{\mathbf{r}}={\mathbf{ϵ}}_{\mathbf{r}}^{\mathrm{IG}}+{\mathbf{ϵ}}_{\mathbf{r}}^{\mathrm{CG}}\sim \mathcal{N}(0,{\mathbf{Q}}_{\mathbf{r}}^{\mathrm{IG}}+{\mathbf{Q}}_{\mathbf{r}}^{\mathrm{CG}})\text{and}$$

$${\mathbf{ϵ}}_{\mathbf{\omega }}={\mathbf{ϵ}}_{\mathbf{\omega }}^{\mathrm{IG}}+{\mathbf{ϵ}}_{\mathbf{\omega }}^{\mathrm{CG}}\sim \mathcal{N}(0,{\mathbf{Q}}_{\mathbf{\omega }}^{\mathrm{IG}}+{\mathbf{Q}}_{\mathbf{\omega }}^{\mathrm{CG}}).$$

$${\mathbf{\Sigma }}_{\mathrm{CG}}={\zeta }^2\left[\begin{array}{cc}{c}_6^2& c_6{c}_7\\ c_6{c}_7& c_7^2\end{array}\right].$$

$$f_{h_g}(h)=$$

Full text

I Introduction

Recently, there has been a growing interest in unmanned aerial vehicles (UAVs) for civil applications, such as delivering wireless access to remote regions or areas where a large number of users is temporarily gathered, e.g., for a football match or a live concert, and permanent infrastructure is not available or is costly to deploy [2, 3]. In particular, UAVs may hover above the desired area and operate as mobile remote radio heads to assist the communication between the users and a central unit (CU) [2].

For these applications, free space optics (FSO) communication has been considered as a promising candidate for fronthauling of the data gathered by the UAVs to the CU [2, 4, 5]. FSO systems offer the large bandwidth needed for data fronthauling, while FSO transceivers are relatively cheap compared to their radio frequency (RF) counterparts and easy to deploy [5, 6]. However, the quality of the FSO link between a hovering UAV and a CU is negatively affected by variations (jitters) of the position and orientation of the UAV, which originate from several sources including dynamic wind load, inherent random air fluctuations in the atmosphere around the UAV, and internal vibrations of the UAV caused by the rotation of its propellers. These variations directly affect the performance of the tracking system, which is responsible for aligning the beam with the receiver lens at the CU [7, 8, 9, 10, 11]. Therefore, one important question is: How well (stable) does a UAV have to maintain its position and orientation in order to achieve a certain FSO link quality? In this paper, we develop a mathematical framework for answering this question by statistically characterizing the geometric and misalignment losses (GML)111 The receiver can only capture the fraction of power that falls on its lens. This phenomenon is known as geometric loss. Moreover, misalignment of the center of the optical beam and the center of the receiver lens further increases the geometric loss. This phenomenon is known as misalignment loss [12]. caused by the random fluctuations of the position and orientation of UAVs.

We note that even for conventional FSO systems with immobile transceivers fixed at building tops, random fluctuations of the positions of the transceivers occur due to building sway, which leads to random GML, known as pointing errors [12, 13, 14]. For this case, corresponding statistical models were derived in [12] and [13]. However, UAV-based FSO systems introduce the following new challenges: i) For conventional FSO links, it is typically assumed that the laser beam is orthogonal with respect to (w.r.t.) the receiver lens plane, as orthogonality maximizes the amount of laser power collected by the photo-detector (PD) located behind the lens [12]. However, orthogonality may not hold for UAV-based FSO communication systems. For example, the position of a UAV may depend on the locations and traffic needs of the users, while the CU may not be able to adjust the orientation of the receiver lens due to limited mechanical capabilities. In addition, using one receiver lens and multiple PDs [15, 16], the CU may receive data from several UAVs having different positions. Hence, it is not possible to orthogonally align the laser beams of all UAVs with the receiver lens plane. ii) Unlike building sway, where the buildings exhibit limited movement due to wind loads and thermal expansion, for UAV-based FSO communication, both the position and orientation of the UAV may fluctuate over time and have to be modeled as random variables (RVs).

UAVs with FSO links have already been considered in the literature [1, 2, 17, 3, 18, 19, 20, 21, 22]. In particular, the authors of [2] discussed the advantages and challenges of FSO fronthauling in UAV-based networks. Moreover, the authors of [17, 3, 18] studied a system consisting of several UAVs that were connected with each other through FSO links. Specifically, the authors of [17] presented a deterministic model for the geometric loss, assuming that the laser beam is always orthogonal to the receiver’s lens plane. To the best of the authors’ knowledge, a statistical model for the GML of a UAV-based FSO channel, which takes into account the fluctuations of the UAV’s position and orientation as well as the non-orthogonality of the laser beam w.r.t. the receiver lens plane, has been reported first in the conference version of this paper [1]. Later on, the authors of [19] derived a statistical model for the GML assuming random UAV positions and orientations, for the special case where the laser beam is orthogonal to the receiver lens plane and the variances of the fluctuations of the position (orientation) are identical for all directions. Moreover, in recent work, the authors of [20, 21, 22] developed statistical channel models for FSO links connecting different UAVs assuming orthogonal beams. In particular, the authors of [20] investigated the outage probability, which was defined as the probability that the transmitting UAV falls out of the receiver’s field-of-view (FoV) due to fluctuations of the orientation of the receiving UAV. Moreover, in [21, 22], adaptive beam control techniques were proposed to cope with the fluctuations of the positions and orientations of the transmitting and receiving UAVs. Unlike [2, 17, 3, 18, 19, 20, 21, 22], we develop statistical channel models for UAV-based FSO communication systems that allow for non-orthogonal laser beams w.r.t. the receiver lens plane and take into account various models for the fluctuations of the position and orientation of the UAV. In particular, this paper makes the following contributions:

• •

We derive the GML for a given position and orientation of the UAV, which we refer to as conditional GML. In particular, since obtaining a closed-form expression for the conditional GML is difficult, if not impossible, we first derive tight lower and upper bounds, and then provide a closed-form approximation based on these bounds.

• •

We consider the following three models for the random fluctuations of the position and orientation of the UAV: i) Independent Gaussian Fluctuations: Since the fluctuations are in general the result of many contributing factors, such as random air fluctuations in the atmosphere around the UAV and internal vibrations of the UAV, by invoking the central limit theorem, we model the resulting fluctuations of the UAV as Gaussian distributed random variables (RVs) [12, 13]. ii) Correlated Gaussian Fluctuations: This more general model allows for correlations that can be the result of e.g. wind causing the UAV to have stronger fluctuations in a certain direction. iii) Correlated Uniform Fluctuations: We also consider uniformly distributed fluctuations which may better model the characteristics of UAV fluctuations with large variance than Gaussian distributed fluctuations [23, 24, 25]. Independent Gaussian, correlated Gaussian, and correlated uniform fluctuations are expected to be suitable models for calm, weakly windy, and strongly windy weather conditions, respectively.

• •

We derive novel statistical models for the GML for each fluctuation scenario. Moreover, we simplify the derived closed-form expressions for some special cases, e.g., when the beam is orthogonal to the receiver lens plane, to obtain further insight.

• •

Based on the developed statistical GML models, we analyze the performance of a UAV-based FSO link in terms of outage probability and ergodic rate. In particular, we assume that the impact of the GML is dominant compared to atmospheric turbulence induced fading. This is a valid assumption when the distance between the UAV and the CU is on the order of several hundred meters, as is validated by simulations in Section VI, cf. Figs. 10 and 10. Next, we derive analytical expressions for the outage probability and ergodic rate of the considered system and analyze their asymptotic behavior for high signal-to-noise ratios (SNRs) for the three statistical GML models.

• •

Simulations are used to validate our derivations and show the impact of the system parameters, e.g., the non-orthogonality of the optical beam, the variance of the fluctuations, and the beam width, on system performance. Our results reveal the existence of a trade-off between the outage probability and the amount of the average (and the maximum) power collected by the receiver. More specifically, when the variance of the fluctuations is large, a wider beam is preferable to avoid outages although this decreases the average (and the maximum) power collected at the receiver. On the other hand, when the variance of the fluctuations is small, a narrower beam is preferable since this increases the amount of power collected by the receiver lens, see Fig. 8.

The remainder of this paper is organized as follows: The system and channel models are presented in Section II. In Section III, we develop the conditional GML model, and in Section IV, we derive statistical GML models for three different fluctuation scenarios. In Section V, we analyze the performance of a UAV-based FSO link using the developed channel models. In Section VI, we present simulation results, and Section VII concludes the paper.

Notations: Boldface lower-case and upper-case letters are reserved for vectors and matrices, respectively. $\mathbbmss{E}\{\cdot\}$, $(\cdot)^{\mathsf{T}}$, and $\|\cdot\|$ denote expectation, the transpose of a matrix, and the $l_{2}$-norm of a vector, respectively. $\mathbb{R}$ and $\mathbb{R}^{+}$ denote the sets of real and positive real numbers, respectively. $\mathbf{I}$ represents the identity matrix and $\mathrm{diag}\{a_{1},\dots,a_{n}\}$ denotes a diagonal matrix with $a_{1},\dots,a_{n}$ on its main diagonal. $\ln(\cdot)$, $\mathrm{erf}(\cdot)$, $Q(\cdot)$, and $Q(\cdot,\cdot)$ denote the natural logarithm, the error function, the Gaussian Q-function, and the first-order Marcum Q-function, respectively. $\mathbf{a}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ is used to indicate that $\mathbf{a}$ is a multivariate Gaussian random vector with mean vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$ and $b\sim\mathcal{U}(a,b)$ means that RV $b$ is uniformly distributed in interval $[a,b]$. Finally, $\mathbf{a}\cdot\mathbf{b}$ and $\mathbf{a}\times\mathbf{b}$ denote the dot and cross products of vectors $\mathbf{a}$ and $\mathbf{b}$, respectively.

II System and Channel Models

II-A System Model

We consider a UAV-based uplink transmission where the UAV communicates with mobile users via an RF multiple-access link (e.g., using sub-6 GHz bands) and with the CU via an FSO fronthaul link, see Fig. 1. The focus of this paper is on the fronthaul communication between the UAV and the CU. In particular, we assume that the UAV is equipped with an aperture FSO transmitter pointing towards the CU, which detects the received optical power. To avoid the requirement of a mechanical adjustment of the orientation of the receiver at the CU, we assume that the receiver is equipped with a lens placed in front of an array of PDs [15, 16, 26], cf. Fig. 1. The lens separates the signals coming from different spatial angles and focuses them on corresponding PDs, respectively. The large overall FoV of the PD array enables the CU to receive data when a moving UAV is at different positions. Moreover, this architecture allows the CU to simultaneously receive data from multiple UAVs which are spatially separated but are still within the receiver FoV [16]. The overall receiver FoV and the FoV of the individual PDs are design parameters which depend on the specific implementation [15, 16]. As mentioned before, the main goal of this paper is to develop a mathematical framework that models the impact of the fluctuations of the position and orientation of a hovering UAV on the FSO channel quality. Therefore, we assume an ideal PD array where one of the PDs is able to collect the entire optical power flux into the receiver lens (albeit with a fixed efficiency/responsivity factor, see Section II-B).

To characterize an object in three dimensions, at most six independent variables are needed: three variables to specify the position of a reference point of the object and another three to quantify its orientation. Next, we define the position and orientation of the UAV and the CU.

II-A1 CU

The CU is a fixed node located at the top of a building222The CU may not be stable due to building sway [12, 13]. Nevertheless, since only the relative movement of UAV and CU affects the FSO channel quality, we assume that the CU is fixed and only the UAV moves. Note that, in practice, the movement of the CU is negligible compared to the movement of the UAV.. Without loss of generality, we choose the center of the receiver lens as the reference point, which is located in the origin of the Cartesian coordinate system $(x,y,z)=(0,0,0)$. This coordinate system is referred to as Coordinate System 1, cf. Fig. 2. Moreover, we assume a circular lens of radius $r_{0}$. Note that it suffices to characterize the plane where the receiver lens lies in order to specify its orientation. Here, without loss of generality, we assume that the lens lies in the $y-z$ plane at $x=0$.

II-A2 UAV

For the communication system under consideration, the parameters that directly affect the FSO channel are the position of the laser source of the UAV and the direction of the laser beam. Therefore, without loss of generality, we refer to them as the position and orientation of the UAV, respectively. Furthermore, we assume that the UAV is in the hovering state. However, in practice, the position and orientation of the UAV are not perfectly constant in the hovering state [27, 28, 29] and thus, they are modeled as RVs. In particular, let $\mathbf{r}=(r_{x},r_{y},r_{z})$ and $\boldsymbol{\omega}=(\theta,\phi,\wp)$ denote the vectors containing the random position and orientation variables of the UAV, respectively. Without loss of generality, in order to simplify the analysis, we define vector $\mathbf{r}$ w.r.t. Coordinate System 1, whereas we use the following coordinate system for $\boldsymbol{\omega}$: For a given vector $\mathbf{r}$, we define Coordinate System 2 with $\mathbf{r}$ as its origin and axes $x^{\prime}$, $y^{\prime}$, and $z^{\prime}$ that are parallel to the $x$, $y$, and $z$ axes of Coordinate System 1, respectively, cf. Fig. 2. We use variables $\theta$ and $\phi$ to determine the direction of the laser beam in a spherical representation of Coordinate System 2. In particular, $\theta\in[0,2\pi]$ denotes the angle between the projection of the beam vector onto the $x^{\prime}-y^{\prime}$ plane and the $x^{\prime}$ axis; and $\phi\in[0,\pi]$ represents the angle between the beam vector and the $z^{\prime}$ axis. The third orientation variable $\wp$ is used to quantify the rotation around the beam vector. This representation of the orientation variables has two advantages. First, variable $\boldsymbol{\omega}$ does not change if position $\mathbf{r}$ changes, i.e., the position and orientation variables are independent. Second, a rotation around the beam line does not affect the signal at the PD assuming rotational beam symmetry. Therefore, the value of $\wp$ is irrelevant for the analysis, and hereafter, for simplicity, we drop $\wp$ and use $\boldsymbol{\omega}=(\theta,\phi)$ as the random vector of the orientation variable.

II-B FSO Channel Model

We assume an intensity modulation/direct detection (IM/DD) FSO system, where the PD responds to changes in the received optical signal power [4]. Moreover, we assume that noise caused by background illumination is the dominant noise source at the PD [4, 30]. Hence, in our model, the noise is independent from the signal. The received signal at the CU is given by

[TABLE]

where $x_{s}\in\mathbb{R}^{+}$ is the transmitted optical symbol (intensity), $n\in\mathbb{R}$ is the zero-mean real-valued additive white Gaussian shot noise with variance $\sigma_{n}^{2}$ caused by background illumination at the CU, and $h\in\mathbb{R}^{+}$ denotes the FSO channel coefficient. Moreover, we assume an average power constraint $\mathbbmss{E}\{x_{s}\}\leq{P}$. The FSO channel coefficient, $h$, is affected by several phenomena and can be modeled as [12]

[TABLE]

where $\eta$ is the responsivity of the PD and $h_{p}$, $h_{a}$, and $h_{g}$ represent the atmospheric loss, atmospheric turbulence induced fading, and GML, respectively. In particular, the atmospheric loss, $h_{p}$, is deterministic and represents the power loss over a propagation path due to absorption and scattering of the light by particles in the atmosphere. It is modeled as [31, 32]

[TABLE]

where $L$ is the distance between the UAV and the CU and $\kappa$ [m*-1*] denotes the weather-dependent attenuation constant of the FSO link. For clear air, haze, light fog, moderate fog, and heavy fog, the typical values of $\kappa$ are $\{0.43,4.2,20,42.2,125\}\times 10^{-3}$ m*-1*, respectively [32].

The atmospheric turbulence, $h_{a}$, is an RV and induced by inhomogeneities in the temperature and the pressure of the atmosphere. It is typically modeled as log-normal (LN) and Gamma-Gamma (GG) distributed RV for weak and moderate-to-strong turbulence conditions [12], respectively. For the considered system, the distance between the UAV and the CU is typically on the order of several hundred meters. In this regime, the atmospheric turbulence is moderate and its impact is negligible compared to that of the GML. To verify this claim, let us consider the pessimistic GG fading model, i.e., $h_{a}\sim\mathcal{GG}(\alpha,\beta)$, with fading parameters $\alpha$ and $\beta$ [33]. In particular, for GG fading, $h_{a}$ is modeled as the product of two independent Gamma random variables $h_{a}^{(1)}\sim\mathcal{G}(\alpha,\alpha)$ and $h_{a}^{(2)}\sim\mathcal{G}(\beta,\beta)$, which represent irradiance fluctuations arising from large- and small-scale turbulences, respectively [33]. Parameters $\alpha$ and $\beta$ are the inverse of the variances of $h_{a}^{(1)}$ and $h_{a}^{(2)}$, respectively, and are given by [13]

[TABLE]

In (II-B), $\sigma_{R}^{2}=1.23C_{n}^{2}k^{7/6}L^{11/6}$ is the Rytov variance, $k=2\pi/\lambda$, where $\lambda$ [m] denotes the optical wavelength, and $C_{n}^{2}\approx C_{0}^{2}\exp\left(-\frac{h_{d}}{100}\right)$ [m${}^{-\frac{2}{3}}$] is the index of refraction structure parameter, where $h_{d}$ is the operating height of the UAV and $C_{0}^{2}=1.7\times 10^{-14}$ m${}^{-\frac{2}{3}}$ is the nominal value of the refractive index at the ground [13]. For typical system parameters, the variance of $h_{a}$, i.e., $\frac{1}{\alpha}+\frac{1}{\beta}+\frac{1}{\alpha\beta}$, is very small (e.g., $3\times 10^{-2}$ for $L=500$ m, $h_{d}=120$ m, and $\lambda=1550$ nm). Therefore, we approximate $h_{a}$ by its mean value, i.e., $h_{a}\approx\mathbbmss{E}\{h_{a}\}=\mathbbmss{E}\{h_{a}^{(1)}\}\mathbbmss{E}\{h_{a}^{(2)}\}=\frac{\alpha}{\alpha}\times\frac{\beta}{\beta}=1$. We verify this assumption by simulations in Section VI, cf. Figs. 10 and 10.

The GML, $h_{g}$, is caused by the divergence of the optical beam between the transmitter and the receiver lens and the misalignment of the laser beam line and the center of the lens [4, 32]. Fluctuations of the position and orientation of the UAV lead to a random GML, $h_{g}$. In the following, we first derive a conditional model for the GML and then, we develop statistical models for independent Gaussian, correlated Gaussian, and uniformly distributed fluctuations.

III The Conditional GML Model

In this section, we derive the channel parameter $h_{g}$ for a given state of the UAV, i.e., for given $\mathbf{r}$ and $\boldsymbol{\omega}$.

III-A Center of the Beam Footprint

The line of the beam can be represented in Cartesian Coordinate System 1 as

[TABLE]

where $\jmath$ is an arbitrary real number and $\mathbf{d}=(d_{x},d_{y},d_{z})$ denotes the beam direction, which can be written as a function of $\theta$ and $\phi$ as

[TABLE]

The center of the beam footprint on the receiver lens can be obtained as the intersection point of the line of the laser beam and the lens plane, $x=0$. Denoting the center of the footprint of the beam on the lens as $\mathbf{b}=(b_{x},b_{y},b_{z})$, then

[TABLE]

III-B Power Density on the Lens Plane

We assume a Gaussian beam, which dictates that the power density distribution across any plane perpendicular to the direction of the wave propagation follows a Gaussian profile [4, 12]. In particular, we consider a perpendicular plane where the distance between the center of the beam footprint on the plane and the laser source is denoted by $L$. Then, the power density for any point on this perpendicular plane with distance $l$ from the center of the beam footprint is given by [12]

[TABLE]

where $w_{L}$ [m] is the beam width at distance $L$ and can be evaluated as

[TABLE]

For the case where the beam propagates in the $x$ direction, $l=\sqrt{\tilde{y}^{2}+\tilde{z}^{2}}$ holds where $\tilde{y}=y-b_{y}$ and $\tilde{z}=z-b_{z}$. In (9), $w_{0}$ [m] denotes the beam waist radius and $\rho(L)=(0.55C_{n}^{2}k^{2}L)^{-3/5}$ [m] is referred to as the coherence length. Recall that for the problem at hand, the plane of the receiver lens is not necessarily orthogonal to the beam direction. For this case, the power density on the lens plane, denoted by $I(y,z)$, is given in the following lemma.

Lemma 1

Under the mild conditions $\|\mathbf{r}\|\gg\|\mathbf{b}\|$ and $\|\mathbf{r}\|\gg\|(y,z)\|$, the power density at point $(y,z)$ on the PD plane is given by

[TABLE]

where $\psi=\sin^{-1}(\sin\phi\cos\theta)$ is the angle between the beam line and the lens plane, $l(\boldsymbol{\omega},y,z)=\sqrt{\rho_{y}\tilde{y}^{2}+\rho_{z}\tilde{z}^{2}+2\rho_{yz}\tilde{y}\tilde{z}}$, $L(\mathbf{r})=\|\mathbf{r}\|$, and $I^{\mathrm{orth}}(\cdot;\cdot)$ is given by (8). Moreover, $\rho_{y}=\cos^{2}\phi+\sin^{2}\phi\cos^{2}\theta$, $\rho_{z}=\sin^{2}\phi$, and $\rho_{yz}=-\cos\phi\sin\phi\sin\theta$.

Proof:

The proof is given in Appendix A. ∎

Note that the conditions under which (10) holds are met in practice since, for typical FSO links, $\|\mathbf{r}\|$ is on the order of several hundred meters, whereas $\|\mathbf{b}\|$ and $\|(y,z)\|$ are on the order of a few centimeters.

III-C GML

The fraction of power collected by the receiver lens, denoted by $h_{g}(\mathbf{r},\boldsymbol{\omega})$, can be obtained by integrating the power density derived in Lemma 1 over the lens area. This leads to

[TABLE]

where $I(y,z)$ is given in (10) and $\mathcal{A}=\left\{(y,z)|y^{2}+z^{2}\leq r_{0}^{2}\right\}$ is the set of $(y,z)$ within the lens area. The exact value of the integral in (11) cannot be obtained in closed form. Instead, in the following theorem, we provide an upper and a lower bound on $h_{g}(\mathbf{r},\boldsymbol{\omega})$.

Theorem 1

Using Lemma 1, $h_{g}(\mathbf{r},\boldsymbol{\omega})$ can be lower and upper bounded by

[TABLE]

*respectively. Here, $u=\|\mathbf{b}\|$ denotes the distance between the origin and the center of the beam footprint, i.e., the misalignment. *

Proof:

The proof is given in Appendix B. ∎

Remark 1

We use Fig. 3 to illustrate the basic idea behind the upper and lower bounds proposed in Theorem 1. In particular, unlike the case where the optical beam is orthogonal to the lens plane and the power density contours are circles [12], when the optical beam is non-orthogonal to the lens plane, the power density contours are rotated ellipses, e.g., the red ellipse in Fig. 3. We have derived the lower bound assuming a footprint that is a rotated ellipse, whose major axis is perpendicular to the line connecting the center of the footprint and the origin, i.e., the green ellipse in Fig. 3. Moreover, for the upper bound, the footprint is a rotated ellipse, whose minor axis is perpendicular to the line connecting the center of the footprint and the origin, i.e., the purple ellipse in Fig. 3. In the special case where the major (minor) axis of the original power density contour is perpendicular to the line connecting the center of the footprint and the origin, the lower (upper) bound is identical to the exact GML.

The integrals in (1) cannot be evaluated in closed-form. Even for the case where the beam line is orthogonal to the lens plane (as is the case for conventional FSO systems [12]), the exact value of $h_{g}(\mathbf{r},\boldsymbol{\omega})$ is cumbersome and provides little insight. Therefore, in [12], the authors proposed an approximation for conventional FSO systems, which was shown to be very accurate for $\frac{w_{L}}{r_{0}}\geq 6$ and has been widely used by other authors subsequently [34, 14, 13, 19]. The proposed bounds in Theorem 1 have two main advantages. First, for the special case where the beam line is orthogonal to the lens plane, the upper and lower bounds coincide with the exact $h_{g}(\mathbf{r},\boldsymbol{\omega})$. Second, the form of the integrals in (1) allows to employ the same technique as in [12, Appendix] in order to obtain accurate approximations. In particular, as shown in detail in Appendix C, we approximate $h_{g}^{\mathrm{low}}(\mathbf{r},\boldsymbol{\omega})$ and $h_{g}^{\mathrm{upp}}(\mathbf{r},\boldsymbol{\omega})$ in (1) with $\widetilde{h}_{g}^{\mathrm{low}}(\mathbf{r},\boldsymbol{\omega})$ and $\widetilde{h}_{g}^{\mathrm{upp}}(\mathbf{r},\boldsymbol{\omega})$, respectively, as follows

[TABLE]

where $t_{1}=\frac{\sqrt{\pi}\mathrm{erf}(\nu_{1})}{2\nu_{1}\exp(-\nu_{1}^{2})}$, $\nu_{1}=\frac{r_{0}}{w_{L}}\sqrt{\frac{\pi}{2}}$, $t_{2}=\frac{\sqrt{\pi}\mathrm{erf}(\nu_{2})}{2\nu_{2}\exp(-\nu_{2}^{2})\sin^{2}\phi\cos^{2}\theta}$, and $\nu_{2}=\nu_{1}|\sin\phi\cos\theta|$. Moreover, $A_{0}$ denotes the maximum fraction of optical power captured by the receiver lens at $u=0$ and is given by

[TABLE]

Note that $A_{0}$ is inversely proportional to $\frac{w_{L}}{r_{0}}$, which means that, as expected, the wider the beam footprint w.r.t. the lens is, the smaller the amount of power that can be collected by the lens. The only difference between the approximated lower and upper bounds in (14a) are the factors $t_{1}$ and $t_{2}$. This motivates us to propose the following approximation for the GML

[TABLE]

where $t\in[t_{1},t_{2}]$. In (17), $h_{g}(\mathbf{r},\boldsymbol{\omega})$ comprises two parts, namely $A_{0}$, which affects the geometric loss and $\exp\left(-\frac{2u^{2}}{tw^{2}_{L}}\right)$, which represents the misalignment attenuation when $u\neq 0$.

In the following, instead of considering the approximate upper and lower bounds in (14a), we employ the general approximation in (17) for statistical analysis. One can choose $t$ in (17) equal to $t_{1}$ and $t_{2}$ to obtain the lower and upper bounds, respectively. Alternatively, $t$ can be chosen as the arithmetic mean $\frac{t_{1}+t_{2}}{2}$ or the geometric mean $\sqrt{t_{1}t_{2}}$ to compromise between the lower and upper bounds. Our results in Section VI show that the approximation in (17) yields an accurate approximation of $h_{g}(\mathbf{r},\boldsymbol{\omega})$ for both $t=\frac{t_{1}+t_{2}}{2}$ and $t=\sqrt{t_{1}t_{2}}$ for the practical range of system parameters, cf. Fig. 4.

IV Statistical Models for the GML

In (17), we provided an approximate closed-form expression for the GML $h_{g}(\mathbf{r},\boldsymbol{\omega})$ for given values of $\mathbf{r}$ and $\boldsymbol{\omega}$. However, in practice, the position and orientation of a hovering UAV fluctuate randomly, and hence, $\mathbf{r}$ and $\boldsymbol{\omega}$ are RVs. In the following, we first present three fluctuation scenarios for RVs $\mathbf{r}$ and $\boldsymbol{\omega}$ and then derive the corresponding statistical GML models.

IV-A Models for the Random Position and Orientation Fluctuations of the UAV

As mentioned above, the position and orientation of the UAV randomly fluctuate over time. In other words, a hovering UAV is not perfectly stable [17]. Therefore, an active control mechanism is needed to persistently keep the laser beam and the receiver lens aligned (see [7, 8, 9, 11]). For ideal tracking, the center of the beam footprint coincides with the center of the receiver lens, i.e., $u=0$. Nevertheless, in practical systems, misalignment due to tracking errors exists for several reasons. For instance, the control system requires some time to compensate for a misalignment or this system is not perfectly accurate, while compensating for a misalignment. Moreover, in UAVs, there is an error associated with wind estimation, i.e., the power and direction of wind[28], and therefore, the impact of wind cannot be fully compensated. In fact, tracking errors exist even in conventional FSO systems, where the transceivers are mounted on top of buildings and misalignment originates from building sway. However, for UAV-based FSO links, such tracking errors are expected to be more severe, due to the inherent instability of hovering UAVs. Therefore, for the development of a channel model for UAV-based FSO links, statistical models for the position and orientation of the UAV are needed.

Let us define vectors $\boldsymbol{\mu}_{\mathbf{r}}=(\mu_{x},\mu_{y},\mu_{z})$ and $\boldsymbol{\mu}_{\boldsymbol{\omega}}=(\mu_{\theta},\mu_{\phi})$, which denote the means of random vectors $\mathbf{r}$ and $\boldsymbol{\omega}$, respectively. Furthermore, we define the zero-mean random vectors $\boldsymbol{\epsilon}_{\mathbf{r}}=({\epsilon}_{x},{\epsilon}_{y},{\epsilon}_{z})$ and $\boldsymbol{\epsilon}_{\boldsymbol{\omega}}=({\epsilon}_{\theta},{\epsilon}_{\phi})$ to model the fluctuations of the position and orientation of the UAV, respectively. Therefore, the position and orientation of the UAV are respectively given by

[TABLE]

Since the GML is a function of $\boldsymbol{\epsilon}_{\mathbf{r}}$ and $\boldsymbol{\epsilon}_{\boldsymbol{\omega}}$, the distribution of $\boldsymbol{\epsilon}_{\mathbf{r}}$ and $\boldsymbol{\epsilon}_{\boldsymbol{\omega}}$ determines the distribution of the GML. Hence, adopting appropriate distributions for fluctuations $\boldsymbol{\epsilon}_{\mathbf{r}}$ and $\boldsymbol{\epsilon}_{\boldsymbol{\omega}}$ is important for developing a realistic statistical model for the GML. Hereby, independent Gaussian, correlated Gaussian, and correlated uniformly distributed fluctuations may serve as models for calm, weakly windy, and strongly windy weather conditions, respectively. In the following, we discuss the first and second moments of RVs $\mathbf{r}$ and $\boldsymbol{\omega}$.

IV-A1 First Moments of RVs

Since the UAV is supposed to hover above the area where the users are located, $\boldsymbol{\mu}_{\mathbf{r}}$ depends on the location of the users as well as the desired operating height of the UAV. Given $\boldsymbol{\mu}_{\mathbf{r}}$, the tracking system of the UAV initially aims to determine $\boldsymbol{\mu}_{\boldsymbol{\omega}}$ such that the beam line intersects with the center of the receiver lens, i.e., the origin, such that the lens collects the maximum power. This leads to

[TABLE]

In other words, $\mathbb{E}\{\mathbf{b}\}=(0,0,0)$.

IV-A2 Second Moments of RVs

The second moments of $\mathbf{r}$ and $\boldsymbol{\omega}$ determine how well the UAV is able to maintain its position and orientation around the mean values $\boldsymbol{\mu}_{\mathbf{r}}$ and $\boldsymbol{\mu}_{\boldsymbol{\omega}}$, respectively. In particular, the smaller the variances of the elements of vectors $\mathbf{r}$ and $\boldsymbol{\omega}$ are, the more stable the UAV is. Hence, we consider the variances of the position and orientation of the UAV as a measure for the stability of the UAV and subsequently evaluate the performance of the FSO fronthaul link in terms of this measure.

IV-B Statistical GML Model for Independent Gaussian Fluctuations

Position and orientation of a hovering UAV fluctuate around their mean values even in calm weather conditions, i.e., in the absence of wind. These fluctuations are the result of many factors such as inherent random air fluctuations in the atmosphere around the UAV and the internal vibrations of the UAV due to e.g. the rotation of its propellers. Hence, invoking the central limit theorem, we model the resulting position and orientation fluctuations of the UAV as Gaussian distributed RVs. Moreover, we assume that the fluctuations are independent. We note that this is inline with the independent Gaussian fluctuations assumed for derivation of the statistical model for the geometric spread and pointing error due to building sway for conventional FSO links [12, 13]. Fluctuations $\boldsymbol{\epsilon}_{\mathbf{r}}$ and $\boldsymbol{\epsilon}_{\boldsymbol{\omega}}$ are modeled as zero-mean Gaussian random vectors, i.e.,

[TABLE]

where the elements of $\boldsymbol{\epsilon}_{\mathbf{r}}^{{\mathrm{IG}}}=({\epsilon}_{x}^{{\mathrm{IG}}},{\epsilon}_{y}^{{\mathrm{IG}}},{\epsilon}_{z}^{{\mathrm{IG}}})$ and $\boldsymbol{\epsilon}^{{\mathrm{IG}}}_{\boldsymbol{\omega}}=({\epsilon}_{\theta}^{{\mathrm{IG}}},{\epsilon}_{\phi}^{{\mathrm{IG}}})$ are independent Gaussian RVs. Moreover, covarinace matrices $\mathbf{Q}_{\mathbf{r}}^{\mathrm{IG}}$ and $\mathbf{Q}_{\boldsymbol{\omega}}^{\mathrm{IG}}$ are defined as $\mathbf{Q}_{\mathbf{r}}^{\mathrm{IG}}=\mathrm{diag}\{\sigma_{x}^{2},\sigma_{y}^{2},\sigma_{z}^{2}\}$ and $\mathbf{Q}_{\boldsymbol{\omega}}^{\mathrm{IG}}=\mathrm{diag}\{\sigma_{\theta}^{2},\sigma_{\phi}^{2}\}$, where $\sigma^{2}_{s},\,\,s\,\in\{x,y,z,\theta,\phi\},$ is the variance of the fluctuation of component $s$.

The PDF of the GML can be derived by combining (17)-(20). Note that in (17), $A_{0}$, $t$, and $u$ are RVs since $A_{0}$ and $t$ depend on RV $\boldsymbol{\omega}$, and $u$ depends on both RVs $\mathbf{r}$ and $\boldsymbol{\omega}$. However, the variances of $A_{0}$ and $t$ are several orders of magnitude smaller than the variance of $u$. The reason for this is that a small variation in $\boldsymbol{\omega}$, e.g., on the order of mrad, has a significant impact on $u=\sqrt{b_{y}^{2}+b_{z}^{2}}$ since the impact of this variation on $b_{y}$ and $b_{z}$ in (7) is scaled by $r_{x}$ which typically has a comparatively large value (on the order of several hundred meters). On the other hand, the impact of variations in $\boldsymbol{\omega}$ on $A_{0}$ and $t$ is not scaled by $r_{x}$. Therefore, the fluctuations of $h_{g}(\mathbf{r},\boldsymbol{\omega})$ are mainly caused by the variations of the misalignment, $u$. Hence, in the following, we assume that the values of $A_{0}$ and $t$ are approximately constant and obtained for the average values of the position and orientation of the UAV, i.e., $\boldsymbol{\mu}_{\mathbf{r}}$ and $\boldsymbol{\mu}_{\boldsymbol{\omega}}$. In Section VI, we confirm this assumption via simulations, cf. Figs. 6-10. In addition, as shown in Appendix D, $u$ follows a Nakagami-q (Hoyt) distribution. Based on (17), the relationship between the PDF of $h_{g}$ and $u$, denoted by $f_{h_{g}}(\cdot)$ and $f_{u}(\cdot)$, respectively, is given by

[TABLE]

In the following theorem, we derive the distribution of $h_{g}$ for small $\sigma^{2}_{s},\,\,s\,\in\{x,y,z,\theta,\phi\}$.

Theorem 2

Assuming $\sigma_{s}\to 0,\,s\in\{x,y,z,\theta,\phi\}$, the PDF of $h_{g}$ is given by

[TABLE]

where $\varpi=\frac{(1+q^{2})tw^{2}_{L}}{4q\Omega}$ is a constant and $I_{0}(\cdot)$ is the zero-order modified Bessel function of the first kind. Moreover, $q=\sqrt{\frac{\min\{\lambda_{1},\lambda_{2}\}}{\max\{\lambda_{1},\lambda_{2}\}}}$ and $\Omega=\lambda_{1}+\lambda_{2}$, where $\lambda_{1}$ and $\lambda_{2}$ are the eigenvalues of matrix $\boldsymbol{\Sigma}_{\mathrm{IG}}$, which is given by

[TABLE]

In (23), $c_{1}=-\tan\mu_{\theta}$, $c_{2}=-\frac{\mu_{x}}{\cos^{2}\mu_{\theta}}$, $c_{3}=\frac{\mu_{x}}{\sin^{2}\mu_{\phi}\cos\mu_{\theta}}$, $c_{4}=-\frac{\mu_{x}\cot\mu_{\phi}\tan\mu_{\theta}}{\cos\mu_{\theta}}$, and $c_{5}=-\frac{\cot\mu_{\phi}}{\cos\mu_{\theta}}$ are constants.

Proof:

The proof is given in Appendix D. ∎

Note that the PDF of $f_{h_{g}}(h)$ in (22) has an indeterminate form at $h=0$. Its value can be found for $q\neq 1$ as

[TABLE]

where for equality $(a)$, we used $\underset{z\to\infty}{\lim}I_{0}(z)=\frac{\exp(z)}{\sqrt{2\pi z}}$ [35, Eq. (9.7.1)]. In fact, (24) shows that for a wider beam and smaller variances of the fluctuations, for which $q\varpi\geq 1$ is met, the channel quality becomes better since the probability of small channel coefficient values approaches zero. On the other hand, having a wide beam reduces the maximum fraction of power collected by the receiver lens, $A_{0}$, cf. (16). Therefore, there is a trade-off between $A_{0}$ and $q\varpi$ when choosing the beam width (beam divergence angle). In the following corollary, we investigate the special case when the beam is perpendicular w.r.t. the lens plane.

Corollary 1

When the laser beam is perpendicular w.r.t. the lens plane, i.e., $\mu_{y}=\mu_{z}=0$, $\mu_{\theta}=\pi$, and $\mu_{\phi}=\pi/2$, $\boldsymbol{\Sigma}_{\mathrm{IG}}$ is given by

[TABLE]

which has eigenvalues $\lambda_{1}=\sigma^{2}_{y}+\mu_{x}^{2}\sigma^{2}_{\theta}$ and $\lambda_{2}=\sigma^{2}_{z}+\mu_{x}^{2}\sigma^{2}_{\phi}$. Hereby, assuming $\sigma^{2}_{y}=\sigma^{2}_{z}\triangleq\sigma^{2}_{p}$ and $\sigma^{2}_{\theta}=\sigma^{2}_{\phi}\triangleq\sigma^{2}_{o}$ leads to $q=1$ and RV $u$ is Rayleigh distributed [12]. Therefore, the PDF of $h_{g}$ simplifies to

[TABLE]

where $\varpi=\frac{tw^{2}_{L}}{4(\sigma^{2}_{p}+\mu^{2}_{x}\sigma^{2}_{o})}$.

Proof:

The proof follows by substitution of $q=1$ and $\Omega=2(\sigma^{2}_{p}+\mu^{2}_{x}\sigma^{2}_{o})$ into (22). ∎

Depending on the value of $\varpi$, the PDF of the GML in (26) shows the following behavior. i) If $\varpi>1$ holds, the probability of small channel coefficients becomes very small, i.e., $\underset{h\to 0}{\lim}\,f_{h_{g}}(h)=0$. As a special case when the UAV is fully stable, i.e., $\sigma_{p}=\sigma_{o}=0$ leading to $\varpi\to\infty$, random fluctuations are not present anymore and the GML becomes a deterministic function of the given position and orientation of the UAV. In other words, the PDF of the GML becomes a Dirac function at $A_{0}$, i.e., $f_{h_{g}}(h)=\delta(h-A_{0})$. ii) If $\varpi=1$ holds, the GML is uniformly distributed in $[0,A_{0}]$, i.e., $f_{h_{g}}(h)=\frac{1}{A_{0}}$. iii) If $\varpi<1$ holds, the channel quality deteriorates and the probability of small channel coefficients becomes very large, i.e., $\underset{h\to 0}{\lim}\,f_{h_{g}}(h)=\infty$. Recall that for the non-orthogonal case in (24), when $q\neq 1$, we have two cases $\underset{h\to 0}{\lim}\,f_{h_{g}}(h)\in\{0,\infty\}$ depending on the value of $q\varpi$; whereas for the orthogonal case, $q=1$, we have three cases $\underset{h\to 0}{\lim}\,f_{h_{g}}(h)\in\{0,\frac{1}{A_{0}},\infty\}$ depending on the value of $\varpi$.

The simplified matrix $\boldsymbol{\Sigma}_{\mathrm{IG}}$ in (25) reveals that the GML is much more sensitive to the variance of the orientation, $\sigma^{2}_{o}$, than to the variance of the position, $\sigma^{2}_{p}$, since $\sigma^{2}_{o}$ is scaled by the average distance between the UAV and the CU, i.e., $\mathbb{E}\{\|\mathbf{r}\|\}=\mu_{x}$. Another interesting observation from (25) is that the variation of the position and orientation of the UAV along the $x$ axis does not affect the GML since $\sigma_{x}^{2}$ does not appear in (25). The reason for this is that, since the beam is orthogonal to the receiver lens plane, the optical beam propagates along the $x$ axis, and therefore, small changes of the position of the UAV in $x$ direction do not affect the power collected by the PD. Finally, we note that despite the presence of both position and orientation fluctuations, (26) has a similar form as the geometric spread and the pointing error in conventional FSO systems where only position fluctuations are present [12, 13].

IV-C Statistical GML Model for Correlated Gaussian Fluctuations

Now, we consider the case where there is a weak wind along a specific direction denoted by $\mathbf{v}=(v_{x},v_{y},v_{z})$. In this scenario, it is expected that the wind causes larger fluctuations of RV $\mathbf{r}$ along the direction of $\mathbf{v}$. Similarly, depending on the geometry of the UAV333For a perfect spherical object, due to symmetry, the force applied on its surface by wind does not create rotational forces. However, for practical UAV geometries, the impact of the wind force will be dominant in a certain direction which causes rotational forces. The exact direction of the rotational force depends on the object geometry and the direction of the wind and cannot be specified a priori., the wind may cause larger fluctuations of $\boldsymbol{\omega}$ in a certain direction, denoted by $\boldsymbol{\tau}=(\tau_{\theta},\tau_{\phi})$. Hence, we model the fluctuations of $\mathbf{r}$ and $\boldsymbol{\omega}$ as correlated Gaussian RVs. Note that the total fluctuations are the result of both independent and correlated Gaussian distributed variations. In particular, the fluctuations are modeled as

[TABLE]

Here, $\boldsymbol{\epsilon}_{\mathbf{r}}^{{\mathrm{CG}}}=(\epsilon_{x}^{{\mathrm{CG}}},\epsilon_{y}^{{\mathrm{CG}}},\epsilon_{z}^{{\mathrm{CG}}})$ denotes a random vector with correlated Gaussian distributed elements that models the fluctuation of the position of the UAV due to the wind. In particular, we model $\boldsymbol{\epsilon}_{\mathbf{r}}^{{\mathrm{CG}}}$ as $\boldsymbol{\epsilon}_{\mathbf{r}}^{{\mathrm{CG}}}=\delta^{\mathrm{G}}\mathbf{v}$, where $\delta^{\mathrm{G}}\sim\mathcal{N}(0,\zeta^{2})$ denotes a zero-mean normal RV with variance $\zeta^{2}$. Moreover, $\mathbf{Q}_{\mathbf{r}}^{{\mathrm{CG}}}=\zeta^{2}\mathbf{v}^{\mathsf{T}}\mathbf{v}$ is the covariance matrix of $\boldsymbol{\epsilon}_{\mathbf{r}}^{{\mathrm{CG}}}$. Similarly, we have $\boldsymbol{\epsilon}_{\boldsymbol{\omega}}^{{\mathrm{CG}}}=(\epsilon_{\theta}^{{\mathrm{CG}}},\epsilon_{\phi}^{{\mathrm{CG}}})=\delta^{\mathrm{G}}\boldsymbol{\tau}$ and its covariance matrix is given by $\mathbf{Q}_{\boldsymbol{\omega}}^{{\mathrm{CG}}}=\zeta^{2}\boldsymbol{\tau}^{\mathsf{T}}\boldsymbol{\tau}$ 444For notational consistency, we assume that the unit of $\delta^{\mathrm{G}}$ is meter. This implies that $\mathbf{v}$ is unitless and the unit of $\boldsymbol{\tau}$ is rad.m*-1*.. We note that similar to the independent Gaussian fluctuation model, for the correlated Gaussian fluctuation model, $u$ follows a Nakagami-q (Hoyt) distribution, cf. Appendix E. In the following theorem, we derive the PDF of $h_{g}$.

Theorem 3

Assuming $\sigma_{s}\to 0,\,s\in\{x,y,z,\theta,\phi\}$ and $\zeta\to 0$, $h_{g}$ follows the PDF in (22) if matrix $\boldsymbol{\Sigma}_{\mathrm{IG}}$ is replaced by $\boldsymbol{\Sigma}_{\mathrm{T}}=\boldsymbol{\Sigma}_{\mathrm{IG}}+\boldsymbol{\Sigma}_{\mathrm{CG}}$, where $\boldsymbol{\Sigma}_{\mathrm{CG}}$ is given by

[TABLE]

Here, $c_{6}=v_{y}+v_{x}c_{1}+\tau_{\theta}c_{2}$, $c_{7}=v_{z}+v_{x}c_{5}+\tau_{\phi}c_{3}+\tau_{\theta}c_{4}$, and constants $c_{1}$-$c_{5}$ are defined in Theorem 2.

Proof:

The proof is given in Appendix E. ∎

In the following, we consider the special case where the impact of the wind on the fluctuations of $\mathbf{r}$ and $\boldsymbol{\omega}$ is dominant, i.e., $\mathrm{Tr}\{\boldsymbol{\Sigma}_{\mathrm{CG}}\}\geq\mathrm{Tr}\{\boldsymbol{\Sigma}_{\mathrm{IG}}\}$. In this case, $u$ follows a one-sided Gaussian distribution given in Appendix F.

Corollary 2

For the special case, where the impact of wind is dominant, the PDF of $h_{g}$ simplifies to

[TABLE]

where $\varpi=\frac{tw^{2}_{L}}{4\zeta^{2}(c_{6}^{2}+c_{7}^{2})}$.

Proof:

The proof is given in Appendix F. ∎

For small channel coefficients, i.e., $h\to 0$, the PDF in (29) has the following transient behavior: If $\varpi\geq 1$, $\underset{h\to 0}{\lim}f_{h_{g}}(h)=0$ and if $\varpi<1$, $\underset{h\to 0}{\lim}f_{h_{g}}(h)=\infty$ which is similar to the behavior in the independent Gaussian case, cf. (24). Moreover, $\underset{h\to A_{0}}{\lim}f_{h_{g}}(h)=\infty$ holds which is different from the independent Gaussian and general correlated Gaussian cases where $\underset{h\to A_{0}}{\lim}f_{h_{g}}(h)=\frac{\varpi}{A_{0}}$ is bounded for $q\neq 0$, cf. (22).

IV-D Statistical GML Model for Correlated Uniform Fluctuations

If strong wind is present, the fluctuations of the position and orientation of the UAV are relatively large compared to those for calm and weakly windy weather conditions. On the other hand, for practical UAVs, it is reasonable to assume that despite being large, the fluctuations are bounded. In this case, assuming a Gaussian distribution for the fluctuations is not appropriate. Instead, the uniform distribution is a better model for the fluctuations of the position and orientation of the UAV. Note that in the absence of prior knowledge, the uniform distribution is a widely-adopted choice for bounded RVs, see e.g., the application of the uniform distribution for robustness analysis in [23], uncertainty analysis in [24], and worst-case analysis in [25]. Similar to Section IV-C, we assume that the UAV position and orientation fluctuations are stronger in a certain direction which may be caused by e.g. wind. Furthermore, we assume that the effect of the wind is the dominant source of the fluctuations. More specifically, the fluctuations of $\mathbf{r}$ and $\boldsymbol{\omega}$ are modeled as

[TABLE]

where $\boldsymbol{\epsilon}_{\mathbf{r}}^{\mathrm{CU}}=(\epsilon_{x}^{{\mathrm{CU}}},\epsilon_{y}^{{\mathrm{CU}}},\epsilon_{z}^{{\mathrm{CU}}})$ denotes the fluctuation of position of the UAV due to strong wind and is given by $\boldsymbol{\epsilon}_{\mathbf{r}}^{\mathrm{CU}}=\delta^{\mathrm{U}}\mathbf{v}$, where $\delta^{\mathrm{U}}\sim\mathcal{U}\left(-\sqrt{3}\xi,\sqrt{3}\xi\right)$ is a zero-mean uniformly distributed RV with variance $\xi$. Similarly, $\boldsymbol{\epsilon}_{\boldsymbol{\omega}}^{\mathrm{CU}}=(\epsilon_{\theta}^{{\mathrm{CU}}},\epsilon_{\phi}^{{\mathrm{CU}}})=\delta^{\mathrm{U}}\boldsymbol{\tau}$ denotes the fluctuation of the orientation of the UAV caused by wind in a specific direction $\boldsymbol{\tau}$. In this case, the misalignment, $u$, follows a uniform distribution, $u\sim\mathcal{U}\left(0,\sqrt{3(c_{6}^{2}+c_{7}^{2})}\,\xi\right)$, cf. Appendix G. The following theorem provides the PDF of the GML for correlated uniform fluctuations.

Theorem 4

Assuming $\xi\to 0$, the PDF of $h_{g}$ is given by

[TABLE]

where $\alpha_{1}=\sqrt{\frac{tw^{2}_{L}}{24(c_{6}^{2}+c_{7}^{2})\xi^{2}}}$ and $h_{1}=A_{0}\exp\left(-\frac{6(c_{6}^{2}+c_{7}^{2})\xi^{2}}{tw^{2}_{L}}\right)$.

Proof:

The proof is given in Appendix G. ∎

We note that PDF $f_{h_{g}}(h)$ assumes large values at $h=A_{0}$, i.e., $\underset{h\to A_{0}}{\lim}f_{h_{g}}(h)=\infty$.

V Performance Analysis

In this section, we analyze the outage probability and ergodic rate of the considered UAV-based FSO system.

V-A Outage Probability

The outage probability is defined as the probability that the SNR, denoted by $\gamma$, falls below a predefined threshold, $\gamma_{\mathrm{thr}}$. For the channel model in (1), the SNR is defined as $\gamma=\eta^{2}h_{p}^{2}h_{g}^{2}\bar{\gamma}$ where $\bar{\gamma}=\frac{P^{2}}{\sigma^{2}_{n}}$ is the transmit SNR. Therefore, the outage probability is obtained as a function of the cumulative distribution function (CDF) of the GML as follows

[TABLE]

where $F_{h_{g}}(\cdot)$ denotes the CDF of the GML. In the following, we derive the outage probability for different fluctuation scenarios.

V-A1 Independent Gaussian Fluctuations

For the case of independent Gaussian fluctuations, using (22), $P_{\mathrm{out}}$ can be written as

[TABLE]

for $0\leq\frac{\sqrt{{\gamma_{\mathrm{thr}}}}}{\eta h_{p}\sqrt{\bar{\gamma}}}\leq A_{0}$, where $a=\frac{1+q}{2q}g$, $b=\frac{1-q}{2q}g$, and $g=\sqrt{\frac{1+q^{2}}{\Omega}\frac{tw^{2}_{L}}{2}\ln\left(\frac{\eta h_{p}A_{0}\sqrt{\bar{\gamma}}}{\sqrt{{\gamma_{\mathrm{thr}}}}}\right)}$ [36]. In the following, we simplify (33) for some special cases.

Corollary 3

For independent Gaussian fluctuations, if the beam is orthogonal to the lens plane, $P_{\mathrm{out}}$ can be obtained as

[TABLE]

Proof:

After integrating the PDF of the GML for this special case (see (26)), the CDF of the GML and hence, $P_{\mathrm{out}}$ can be obtained as in (34). ∎

Eq. (34) reveals that the diversity gain of the FSO link is

[TABLE]

Corollary 4

For high SNRs, i.e., for large values of arguments $a$ and $b$, (33) can be simplified as [37]

[TABLE]

where $a_{t}=\frac{2\sqrt{2\Omega}b_{t}^{\frac{(1+q^{2})tw_{L}^{2}}{4\Omega}}}{\sqrt{\pi(1-q^{4})tw_{L}^{2}}}$ and $b_{t}=\frac{\sqrt{{\gamma_{\mathrm{thr}}}}}{\eta h_{p}A_{0}}$.

Proof:

The proof is given in Appendix H. ∎

Based on (36), the diversity gain is

[TABLE]

since, as $\bar{\gamma}\to\infty$, the impact of the logarithmic term $\left[\ln\left(\frac{\bar{\gamma}}{b_{t}^{2}}\right)\right]^{-\frac{1}{2}}$ in (36) becomes negligible compared to that of the polynomial term $\left(\frac{1}{\bar{\gamma}}\right)^{\frac{(1+q^{2})tw_{L}^{2}}{8\Omega}}$.

V-A2 Correlated Gaussian Fluctuations

In this case, if both independent and correlated Gaussian fluctuations are present, then the same expression for the outage probability holds, as for the independent Gaussian scenario, cf. (33)-(37). If the impact of wind on the fluctuations is dominant, we obtain

[TABLE]

Based on (38) and using the same approximation for the Gaussian Q-function as in Appendix H, the diversity gain can be obtained as

[TABLE]

V-A3 Correlated Uniform Fluctuations

For uniform distributed fluctuations,

[TABLE]

Note that, in this case, the outage probability is zero if the transmit SNR is larger than a critical value, $\bar{\gamma}_{\mathrm{crt}}$, i.e., $\bar{\gamma}\geq\bar{\gamma}_{\mathrm{crt}}$, where

[TABLE]

V-B Ergodic Rate

For an IM/DD FSO channel, the capacity is not known. Nevertheless, in [38, Eq. 26], the following ergodic rate has been shown to be achievable

[TABLE]

where $c=\frac{e}{2\pi}\eta^{2}h_{p}^{2}\bar{\gamma}$. In the following, we analyze the ergodic rate at high SNR. In particular, for high SNR, we have

[TABLE]

where $\bar{R}_{\max}$ is the maximum achievable ergodic rate without misalignment, i.e., $u=0$, $\Delta\bar{R}_{g}$ is the loss in ergodic rate due misalignment, and $\mathbb{E}_{u}\{u^{2}\}$ denotes the expected value of the squared misalignment, i.e., $u^{2}$. Note that $\Delta\bar{R}_{g}$ depends on the distribution of the fluctuations but $\bar{R}_{\max}$ is independent of it and only depends on mean value of the UAV’s position and orientation, the beam width (beam divergence angle), the transmit SNR, as well as the area and responsivity of the PD. In the following, we evaluate $\Delta\bar{R}_{g}$ for the considered independent/correlated Gaussian and correlated uniform fluctuation models.

Corollary 5

For the considered fluctuation models, $\Delta\bar{R}_{g}$ in bits/symbol is given by

[TABLE]

Proof:

For both independent and correlated Gaussian fluctuation models, $u$ follows a Hoyt distribution. Hence, its second moment is $\mathbb{E}_{u}\{u^{2}\}=\Omega$ [39, Eq. (2.12)] with $\Omega=\lambda_{1}+\lambda_{2}$. For the correlated Gaussian case when the effect of wind is dominant, i.e., $\mathrm{Tr}\{\boldsymbol{\Sigma}_{\mathrm{CG}}\}\geq\mathrm{Tr}\{\boldsymbol{\Sigma}_{\mathrm{IG}}\}$, $\Omega=\lambda_{1}=(c_{6}^{2}+c_{7}^{2})\zeta^{2}$, since $\lambda_{2}=0$. On the other hand, for the uniform fluctuation model, $u$ follows a uniform distribution with second-order moment $\mathbb{E}_{u}\{u^{2}\}=(c_{6}^{2}+c_{7}^{2})\xi^{2}$. Substituting $\mathbb{E}_{u}\{u^{2}\}$ into $\Delta\bar{R}_{g}$ in (43) leads to (44) and concludes the proof. ∎

As can be observed from (44), the rate loss due to misalignment, $\Delta\bar{R}_{g}$, depends on the stability of the UAV through variables $\lambda_{1}$, $\lambda_{2}$, $\zeta$ or $\xi$. Thereby, $\Delta\bar{R}_{g}$ decreases and, as a result, $\bar{R}$ increases as the UAV becomes more stable.

VI Simulations and Discussions

In order to quantify the non-orthogonality of the beam w.r.t. the lens plane, we express the mean position of the UAV, $\boldsymbol{\mu}_{\mathbf{r}}$, in spherical coordinates as $(L,\alpha_{d},\beta_{d})$, i.e., $r_{x}=L\sin\beta_{d}\cos\alpha_{d}$, $r_{y}=L\sin\beta_{d}\sin\alpha_{d}$, and $r_{z}=L\cos\beta_{d}$. Recall that for a given $\boldsymbol{\mu}_{\mathbf{r}}$, the $\boldsymbol{\mu}_{\boldsymbol{\omega}}$ can be obtained from (IV-A1). Unless stated otherwise, the default values of the parameters used for the simulations are: $(\alpha_{d},\beta_{d})=(\frac{\pi}{8},\frac{5\pi}{8})$, $L=500$ m, $h_{d}=120$ m, $\lambda=1550$ nm, $r_{0}=10$ cm, $w_{L}=30$ cm, $\mathbf{v}=\frac{(3,1,2)}{\|(3,1,2)\|}=(0.8,0.27,0.53)$, and $\boldsymbol{\tau}=\frac{1}{L}\frac{(1,2)}{\|(1,2)\|}=\frac{1}{L}(0.44,0.9)$ [12, 13]. Moreover, the simulation results reported in Figs. 6-10 were obtained based on Monte Carlo simulations and $10^{6}$ realizations of RVs $\mathbf{r}$ and $\boldsymbol{\omega}$. These simulation results are used to verify the accuracy of the assumptions made throughout the paper.

First, we study the impact of non-orthogonality of the beam on the conditional GML and investigate the accuracy of the bounds proposed in Theorem 1, their corresponding approximations in (14a), and the proposed approximation for the GML in (17). For Fig. 4, the UAV is located in the $x-y$ plane ($\beta_{d}=\pi/2$) at distance $L$ from the receiver lens and its position on the perimeter of a semicircle with radius $L$ is varied via angle $\alpha_{d}$ w.r.t. the $x$ axis. In particular, in this figure, we show the conditional GML $h_{g}$ vs. $\alpha_{d}$ for two cases of misalignment vector $\mathbf{b}$, namely i) $\mathbf{b}=(0,0,0)$ cm which implies that there is no misalignment, i.e., $u=0$ cm, and ii) $\mathbf{b}=(0,10,10)$ cm which implies that the center of the beam footprint is outside the receiver lens, i.e., $u\approx 14$ cm$>r_{0}$. The curve for $u=0$ cm shows the maximum fraction of power, $A_{0}$, that is collected by the receiver lens for different values of $\alpha_{d}$, and therefore, the bounds in (14a) and the expression proposed for the conditional GML in (17) become identical. Non-zero misalignment causes an additional attenuation \exp\big{(}\frac{-2u^{2}}{tw_{L}^{2}}\big{)} to $A_{0}$, cf. (17). In Fig. 4, this attenuation caused by $14$ cm misalignment is the gap between the curves for $u=0$ and $u=14$ cm at any given $\alpha_{d}$. Moreover, the differences between the maximum value of each curve and any other point on that curve show the loss caused by non-orthogonality of the beam w.r.t. the lens plane. Furthermore, as the beam becomes more non-orthogonal w.r.t. the lens plane, i.e., as $|\alpha_{d}|$ increases, the channel coefficient $h_{g}$ decreases and approaches zero when the beam is parallel to the lens, i.e., $|\alpha_{d}|=\pi/2$. We note that for the practical operating regime of the hovering UAV, $|\alpha_{d}|\leq\pi/4$ holds. Interestingly, in this regime, the proposed approximations for $h_{g}$ in (17), using either the arithmetic or geometric means for $t$, are in good agreement with the simulation results. Besides, in this regime, the loss due to non-orthogonality is small (e.g., $3\times 10^{-2}\approx-1.5$ dB for $u=14$ cm).

Next, in Figs. 6-8, we study the effect of the random fluctuations of the position and orientation of the UAV on the GML and investigate the accuracy of the statistical models developed for different fluctuation scenarios. Specifically, in Fig. 6, the PDF of the GML is plotted assuming independent Gaussian fluctuations of the position and orientation of the UAV with standard deviations (SDs) $(\sigma_{x},\sigma_{y},\sigma_{z})=\sigma r_{0}(0.8,0.27,0.53)$ and $(\sigma_{\theta},\sigma_{\phi})=\frac{\sigma r_{0}}{L}(0.44,0.9)$. Thereby, $\sigma$ controls the SDs and $\frac{1}{L}$ normalizes the orientation fluctuations w.r.t. the distance. In this figure, we plot the PDF of $h_{g}$ for $\sigma=0.5$ and $\sigma=1$ and orthogonal ($(\alpha_{d},\beta_{d})=(0,\frac{\pi}{2})$) and non-orthogonal ($(\alpha_{d},\beta_{d})=(\frac{\pi}{8},\frac{5\pi}{8})$) beams w.r.t. the lens plane. As can be observed from Fig. 6, the analytical statistical model proposed in (22) is in perfect agreement with the histogram obtained based on (11). This agreement also validates our assumption in Section IV.B that the main cause of randomness in the GML is the misalignment $u$ and variables $A_{0}$ and $t$ are practically constant compared to $u$. Moreover, the PDFs for the orthogonal beam, i.e., $(\alpha_{d},\beta_{d})=(0,\frac{\pi}{2})$, assume non-zero values at larger $h_{g}(\mathbf{r},\boldsymbol{\omega})$ compared to those for the non-orthogonal beam, i.e., $(\alpha_{d},\beta_{d})=(\frac{\pi}{8},\frac{5\pi}{8})$, since an additional attenuation is caused by the non-orthogonality of the beam. In addition, for larger $\sigma$, the UAV becomes less stable and the probability of smaller channel coefficients increases. Hence, the corresponding PDFs become more heavy tailed.

Fig. 6 shows the PDF of the GML for the special case considered in Corollary 1 and (26), where the beam is orthogonal to the lens plane and $\sigma_{x}=\sigma_{y}=\sigma_{z}=\sigma_{p}$ and $\sigma_{\theta}=\sigma_{\phi}=\sigma_{o}$ hold, for different variances of the position and orientation fluctuations, i.e., $\sigma_{p}=L\sigma_{o}=r_{0}(\frac{3}{4},1,2)$. It is observed from Fig. 6 that by increasing $\sigma_{p}$ and $\sigma_{o}$, the probability of smaller channel coefficients increases and the PDF becomes more heavy tailed. More specifically, by increasing $\sigma_{p}$ and $\sigma_{o}$, the value of $\varpi$ decreases (see the values of $\varpi$ for different $\sigma_{p}$ and $\sigma_{o}$ in the figure) and for $\varpi<1$, $\lim_{h\to 0}f_{h_{g}}(h)=\infty$ occurs which is consistent with our analytical results in (26).

In Fig. 8, we investigate the case of correlated Gaussian fluctuations. To study the impact of correlation, we use identical variances for the positions (orientations), i.e., the main diagonal entries of $\mathbf{Q}_{\mathbf{r}}^{\mathrm{IG}}$ and $\mathbf{Q}_{\mathbf{r}}^{\mathrm{CG}}$ ($\mathbf{Q}_{\boldsymbol{\omega}}^{\mathrm{IG}}$ and $\mathbf{Q}_{\boldsymbol{\omega}}^{\mathrm{CG}}$), for the independent and correlated Gaussian scenarios are identical. In other words, we set $(\sigma_{x}^{2},\sigma_{y}^{2},\sigma_{z}^{2})=\zeta^{2}(v_{x}^{2},v_{y}^{2},v_{z}^{2})$ and $(\sigma_{\theta}^{2},\sigma_{\phi}^{2})=\zeta^{2}(\tau_{\theta}^{2},\tau_{\phi}^{2})$. In Fig. 8, we plot the PDF of the GML for $\zeta=2r_{0}$, $\mathbf{v}=\frac{(3,4,5)}{\|(3,4,5)\|}=(0.42,0.56,0.7)$, $\boldsymbol{\tau}=(0,0)$, and a non-orthogonal beam w.r.t. the lens plane. It is observed from this figure that for correlated Gaussian fluctuations (Crl. Gauss.), where only the effect of wind is considered, the probabilities of both small and large values for channel coefficient $h_{g}$ are higher compared to the case when the fluctuations are independent (Ind. Gauss.). In particular, the PDF assumes large values at $A_{0}$, which is expected based on our analytical results, cf. (29), and is also large for small values of $h_{g}$ since $\varpi<1$ holds for the set of parameters adopted for this figure (see the value of $\varpi$ in the figure). Furthermore, combined independent and correlated fluctuations (Ind. & Crl. Gauss.) cause the PDF to have larger values for smaller channel coefficients, i.e., it becomes more heavy tailed. Particularly, for the set of parameters adopted in this figure, the PDF for combined independent and correlated Gaussian fluctuations becomes very large at values close to zero since $q\varpi<1$ holds (see the value of $q\varpi$ in the figure), cf. (24).

In Fig. 8, we show the CDF of the GML for correlated uniformly distributed fluctuations for different beam widths, $w_{L}\in\{3r_{0},4r_{0}\}$, and $\xi$, $\xi\in\{3r_{0},4r_{0}\}$. We observe that the simulation and analytical results are in perfect agreement. Moreover, unlike for Gaussian fluctuations, for uniformly distributed fluctuations, the probability of channel coefficients that are smaller than a certain value, $h_{g}\leq h_{1}$, is zero, which is expected based on the analytical expression for the PDF in (31). Comparing the curve for $\xi=3r_{0}$ and the respective curve for $\xi=4r_{0}$ shows that for the larger $\xi$, since the UAV becomes less stable, the channel quality deteriorates, i.e., the value of $h_{1}$ for the CDF for $\xi=4r_{0}$ is smaller than that for $\xi=3r_{0}$. Comparing the curves for different beam widths and a given $\xi$ reveals that, for the wider beam, $w_{L}=4r_{0}$, the maximum fraction of power that is collected at the PD, $A_{0}$, is smaller than $A_{0}$ for $w_{L}=3r_{0}$. On the other hand, given a threshold, the wider the beam is, the smaller the outage probability becomes. For instance, for $\xi=4r_{0}$ in Fig. 8, assuming that an outage occurs when $h_{g}\leq 0.03$, we have $A_{0}=0.16$ and $P_{\mathrm{out}}=0.1$ for $w_{L}=3r_{0}$ and $A_{0}=0.1$ and $P_{\mathrm{out}}=0.02$ for $w_{L}=4r_{0}$. This observation illustrates the trade-off between $A_{0}$ and $P_{\mathrm{out}}$ for different beam widths.

In Figs. 10 and 10, we study the performance of a single UAV-based FSO link in terms of its outage probability and ergodic rate for clear weather conditions, i.e., $\kappa=0.43\times 10^{-3}$ m*-1*. In particular, for Figs. 10 and 10, we assume a non-orthogonal beam w.r.t. the lens plane and independent Gaussian, correlated Gaussian, and correlated uniformly distributed fluctuations for the position and orientation of the UAV. For the independent Gaussian fluctuations, we adopt the same covariance matrices for the position and orientation fluctuations as for Fig. 8. Moreover, we set $\zeta=\xi$ to $r_{0}$ and $\gamma_{\mathrm{thr}}$ is given by $\gamma_{\mathrm{thr}}=\frac{2\pi}{e}2^{2R_{\mathrm{thr}}-1}$, where $R_{\mathrm{thr}}$ is the transmission rate. In Fig. 10, we depict the outage probability vs. SNR ($\bar{\gamma}$) assuming $R_{\mathrm{thr}}=0.5$ bit/symbol in the presence and absence of GG distributed turbulence. We observe that simulation results and analytical results are in perfect agreement for all considered fluctuation models. Furthermore, the gap between the curves with and without GG turbulence is negligible for Gaussian fluctuations which was expected for link lengths, $L$, on the order of several hundred meters, see also the zoomed out part of the figure. The gap is also small for uniform fluctuations for small and medium SNRs but it becomes larger for high SNRs since in this case, for $\bar{\gamma}>\bar{\gamma}_{\mathrm{crt}}$, GG turbulence is the only fading left and the GML is not present anymore, cf. (40). Finally, Fig. 10 confirms the accuracy of the asymptotic outage expression given in (36) for Gaussian fluctuations at high SNRs.

In Fig. 10, we plot the ergodic rate vs. SNR ($\bar{\gamma}$) for the same fluctuation models as considered in Fig. 10. First, the difference between the ergodic rates obtained from simulation with and without GG fading is very small, see also the zoomed out part of the figure. This confirms that the impact of GG turbulence on UAV-based FSO links with lengths on the order of several hundred meters is negligible. Moreover, it is observed that the simulated ergodic rates approach the analytical asymptotic ergodic rates at high SNR for all three fluctuation scenarios. Furthermore, since the asymptotic ergodic rates for different fluctuation models differ only in their constant rate losses $\Delta\bar{R}_{g}$, cf. (43), they approach one another at high SNR. For instance, for the considered set of parameters, the rate losses for independent Gaussian, correlated Gaussian, and correlated uniform fluctuations are given by $\Delta\bar{R}_{g}\in\{0.83,0.9,0.9\}$ bits/symbol, respectively. Moreover, since the $\Delta\bar{R}_{g}$ for the correlated Gaussian and correlated uniform scenarios are identical (since $\zeta=\xi$, cf. (44)), the respective ergodic rates are equal.

VII Conclusions

In this paper, we derived novel statistical models for the FSO fronthaul channel of UAV-based communication systems, by taking into account the non-orthogonality of the laser beam and the random fluctuations of the position and orientation of the UAV. We first modeled the GML as a function of a given position and orientation of the UAV and derived a conditional model. Next, we developed statistical models for the GML assuming independent Gaussian, correlated Gaussian, and correlated uniformly distributed fluctuations for the position and orientation of the UAV which may reflect calm, weakly windy, and strongly windy weather conditions, respectively. Based on the aforementioned channel models, we further analyzed the performance of the UAV-based FSO link in terms of its outage probability and ergodic rate and derived corresponding asymptotic expressions for the high SNR regime. Simulation results validated the presented analysis and revealed important insights for system design. For example, for correlated fluctuations, the probabilities of both small and large values of the channel coefficient are larger than for independent fluctuations. This characteristic leads to a higher outage probability for correlated fluctuations compared to independent fluctuations. Although the specific requirements on the stability of the UAV depend on the targeted application, our simulation results suggested that for high quality FSO links, the standard deviations of the residual fluctuations of the position and orientation of the UAV should not exceed a few $r_{0}$ ($\approx$ tens of cm) and a few $r_{0}/L$ ($\approx$ below one mrad), respectively. We note that the requirements on the permitable fluctuations can be relaxed at the cost of higher power consumption if a wider beam is employed. In particular, when the variance of the fluctuations is relatively large, e.g., due to windy weather conditions or a UAV with low stability, a wider beam is preferable to avoid outages although this causes the average (and the maximum) collected power to decrease. On the other hand, when the variance of the fluctuations is small, a narrower beam is preferable since it increases the amount of power collected by the PD.

Appendix A

$I(y,z)\mathrm{d}y\mathrm{d}z$ is the fraction of power collected in the infinitesimally small area $\mathrm{d}y\mathrm{d}z$, i.e., $\mathrm{d}y\to 0$ and $\mathrm{d}z\to 0$, around the point $(0,y,z)$. Moreover, we use the fact that any point $(0,y,z)$ in the lens plane is also located in another plane which is perpendicular to the beam line. Therefore, power $I(y,z)\mathrm{d}y\mathrm{d}z$ can be obtained as $I(y,z)\mathrm{d}y\mathrm{d}z=I^{\mathrm{orth}}(l;L)\sin\psi\mathrm{d}y\mathrm{d}z$, where $I^{\mathrm{orth}}(l;L)$ is given in (8) and $\psi$ is the angle between the beam line and the lens plane which is found from the inner product of the beam direction and a vector orthonormal to the lens plane, i.e., $(1,0,0)$ in the $x$ direction, as follows

[TABLE]

Here, $\mathbf{d}$ is the direction of the beam given in (6) and we exploited $\|\mathbf{d}\|=1$. Next, we find distances $L$ and $l$. In fact, $l$ is the distance between point $(0,y,z)$ and the beam line in (6). In general, the distance between a point, $\mathbf{p}$, and a line specified by direction vector $\mathbf{u}$ and a given point, $\mathbf{q}$, on the line can be obtained as $l=\frac{\|(\mathbf{p}-\mathbf{q})\times\mathbf{u}\|}{\|\mathbf{u}\|}$. For the problem at hand, we choose $\mathbf{p}=(0,y,z)$, $\mathbf{u}=\mathbf{d}$, and $\mathbf{q}=\mathbf{b}$, which leads to

[TABLE]

where we exploited $\|\mathbf{d}\|=1$, replaced $\mathbf{d}$ with (5), introduced $\tilde{y}=y-b_{y}$ and $\tilde{z}=z-b_{z}$, and used $\rho_{y}$, $\rho_{z}$, and $\rho_{yz}$ given in Lemma 1. Moreover, $L$, the distance between the perpendicular plane w.r.t. the laser beam that contains point $(0,y,z)$ and the laser source can be bounded as

[TABLE]

where the extreme cases occur if the beam line is parallel to the $y-z$ plane. In particular, we can safely assume that $\|\mathbf{r}-\mathbf{b}\|\pm\sqrt{\tilde{y}^{2}+\tilde{z}^{2}}\approx\|\mathbf{r}\|$ holds since the distance between the UAV and the CU, i.e., $\|\mathbf{r}\|$, is much larger than $\|\mathbf{b}\|$ and $\sqrt{\tilde{y}^{2}+\tilde{z}^{2}}$. Therefore, by substituting $L\approx\|\mathbf{r}\|$ and (46) into (8) and using $I(y,z)\mathrm{d}y\mathrm{d}z=I^{\mathrm{orth}}(l;L)\sin\psi\mathrm{d}y\mathrm{d}z$ and (45), we obtain (10) which completes the proof.

Appendix B

In the $y-z$ plane, the contours of power density, $I(y,z)=\bar{I}$, form ellipsoids given by

[TABLE]

These ellipsoids are centered at point $(b_{y},b_{z})$ and rotated by angle $\frac{1}{2}\tan^{-1}$ \big{(}\frac{2\rho_{yz}}{\rho_{y}-\rho_{z}}\big{)} counterclockwise. They have minor and major axis lengths of $2\sqrt{{\rho_{\min}}{d}}$ and $2\sqrt{{\rho_{\max}}{d}}$, respectively, where

[TABLE]

$\rho_{\min}$ and $\rho_{\max}$ can be further simplified using the definition of $\rho_{y}$, $\rho_{z}$, and $\rho_{yz}$ in Lemma 1 as $\rho_{\min}=1$ and $\rho_{\max}=\frac{1}{\sin^{2}\phi\cos^{2}\theta}$.

In order to obtain the lower and upper bounds for $h_{g}(\mathbf{r},\boldsymbol{\omega})$ in Theorem 1, we substitute the contour in (48) by two rotated elliptic contours which have the same axis lengths $\rho_{\min}=1$ and $\rho_{\max}=\frac{1}{\sin^{2}\phi\cos^{2}\theta}$; however, their main axes are either perpendicular or parallel to the line connecting $(b_{y},b_{z})$ and the origin, respectively, see Fig. 3. Moreover, without loss of generality, we can define a new coordinate system by rotating the $y$ and $z$ axes by angle $\tan^{-1}(\frac{b_{z}}{b_{y}})$ such that the center of the ellipsoid in (48) lies on the rotated $y$ axis, i.e., the center becomes $(u,0)$ in the new coordinate system, where $u=\sqrt{b_{y}^{2}+b_{z}^{2}}$. Note that the circular lens has the same description in the new and the old coordinate systems. This leads to lower and upper bounds $h_{g}^{\mathrm{low}}(\mathbf{r},\boldsymbol{\omega})$ and $h_{g}^{\mathrm{upp}}(\mathbf{r},\boldsymbol{\omega})$, respectively, as given in Theorem 1 and completes the proof.

Appendix C

The following integral was approximated in [12, Appendix]

[TABLE]

where $A_{0}=[\mathrm{erf}(\nu)]^{2}$, $t=\frac{\sqrt{\pi}\mathrm{erf}(\nu)}{\sqrt{2}\nu\exp(-\nu^{2})}$, and $\nu=\frac{\sqrt{\pi}r_{0}}{\sqrt{2}w_{L}}$. In (C), equality $(a)$ follows from approximating the circular PD, i.e., $(y,z)\in{\mathcal{A}}$, by a square lens of equal area, i.e., (y,z)\in\bar{\mathcal{A}}\triangleq\Big{\{}(y,z)\mid y,z\in\big{[}-\frac{\sqrt{\pi}r_{0}}{2} ,\frac{\sqrt{\pi}r_{0}}{2}\big{]}\Big{\}} and equality $(b)$ is obtained using the Taylor series of function $\exp(\cdot)$. In the following, we use (C) to approximate $h_{g}^{\mathrm{low}}$ and $h_{g}^{\mathrm{upp}}$ in (1) with $\widetilde{h}_{g}^{\mathrm{low}}$ and $\widetilde{h}_{g}^{\mathrm{upp}}$, respectively. We derive $\widetilde{h}_{g}^{\mathrm{low}}$ since obtaining $\widetilde{h}_{g}^{\mathrm{upp}}$ follows similar steps. By approximating $\mathcal{A}$ with $\bar{\mathcal{A}}$ and defining new variable $\hat{z}=\sin\phi\cos\theta z$ for $h_{g}^{\mathrm{low}}$ in (1a), $\widetilde{h}_{g}^{\mathrm{low}}$ is obtained as

[TABLE]

where \hat{\mathcal{A}}\triangleq\Big{\{}({y},\hat{z})\mid{y}\in\big{[}-\frac{\sqrt{\pi}r_{0}}{2},\frac{\sqrt{\pi}r_{0}}{2}\big{]},\hat{z}\in \big{[}-\frac{\sqrt{\pi}|\sin\phi\cos\theta|r_{0}}{2},\frac{\sqrt{\pi}|\sin\phi\cos\theta|r_{0}}{2}\big{]}\Big{\}}. The integral in (C) is similar to the second integral in (C) except that $\hat{\mathcal{A}}$ corresponds to a rectangular area whereas $\bar{\mathcal{A}}$ is a square area. Using a similar technique as the one used in [12, Appendix], we approximate (C) as in (14a). This completes the proof.

Appendix D

In the following, we first determine the PDF of $u$ and subsequently obtain the PDF of $h_{g}$ from (21). To do so, we first simplify the expressions for $b_{y}$ and $b_{z}$ in (7) by replacing $\tan\theta$, $\cot\phi$, and $\frac{1}{\cos\theta}$ by their respective Taylor series (assuming $\epsilon_{\theta}=\theta-\mu_{\theta}$ and $\epsilon_{\phi}=\phi-\mu_{\phi}$, cf. (18), are very small) and then relating the PDF of $u$ to that of $b_{y}$ and $b_{z}$ exploiting $u=\sqrt{b_{y}^{2}+b_{z}^{2}}$. In particular, we obtain

[TABLE]

where constants $c_{1}$-$c_{5}$ are given in Theorem 2. To obtain (D), we drop the terms with orders higher than one, e.g., $\epsilon_{\theta}\epsilon_{\phi}$. We note that (D) is valid for all considered fluctuation models. Now, assuming independent Gaussian fluctuations, we add superscript IG to $\epsilon_{s},\,\,s\in\{x,y,z,\theta,\phi\}$ and $(b_{y},b_{z})$ in (D), cf. (20). Since $b_{y}^{\mathrm{IG}}$ and $b_{z}^{\mathrm{IG}}$ are sums of Gaussian RVs, they are Gaussian distributed, too. However, $b_{y}^{\mathrm{IG}}$ and $b_{z}^{\mathrm{IG}}$ are correlated since $\epsilon_{x}^{\mathrm{IG}}$ and $\epsilon_{\theta}^{\mathrm{IG}}$ appear in the expressions for both. The joint distribution of $b_{y}^{\mathrm{IG}}$ and $b_{z}^{\mathrm{IG}}$ is a bivariate Gaussian distribution $(b_{y}^{\mathrm{IG}},b_{z}^{\mathrm{IG}})\triangleq\mathbf{b}_{yz}^{\mathrm{IG}}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{\mathrm{IG}})$ where $\boldsymbol{\Sigma}_{\mathrm{IG}}$ is given in (23). Let $\boldsymbol{\Sigma}_{\mathrm{IG}}=\mathbf{U}\boldsymbol{\Lambda}\mathbf{U}^{\mathsf{T}}$ be the eigenvalue decomposition of $\boldsymbol{\Sigma}_{\mathrm{IG}}$ where $\boldsymbol{\Lambda}$ is a diagonal matrix with elements $\lambda_{1}$ and $\lambda_{2}$ and $\mathbf{U}$ is a unitary matrix, i.e., $\mathbf{U}^{\mathsf{T}}\mathbf{U}=\mathbf{I}$. Using these definitions, it is easy to show that $\mathbf{b}_{yz}^{\mathrm{IG}}\sim\mathbf{g}\mathbf{U}^{\mathsf{T}}$ where $\mathbf{g}=(g_{y},g_{z})\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Lambda})$. Now, we can express $u$ in terms of $\mathbf{g}$ as follows

[TABLE]

Since $g_{y}$ and $g_{z}$ are independent zero-mean Gaussian RVs with non-identical variances, $u$ follows a Hoyt (Nakagami-q) distribution with PDF $f_{u}(u)=\frac{1+q^{2}}{q\Omega}u\exp\left(-\frac{(1+q^{2})^{2}}{4q^{2}\Omega}u^{2}\right)I_{0}\left(\frac{1-q^{4}}{4q^{2}\Omega}u^{2}\right)$, where $q=\sqrt{\frac{\min\{\lambda_{1},\lambda_{2}\}}{\max\{\lambda_{1},\lambda_{2}\}}}$ and $\Omega=\lambda_{1}+\lambda_{2}$ [36, 13]. Substituting $f_{u}(u)$ into (21), the PDF of the GML can be obtained as in (22). This completes the proof.

Appendix E

For correlated Gaussian fluctuations, we replace $\epsilon_{s},\,\,s\in\{x,y,z,\theta,\phi\}$ in (D) with $\epsilon_{s}^{\mathrm{IG}}+\epsilon_{s}^{\mathrm{CG}}$ according to the definition in (IV-C). After eliminating the terms of order higher than one, similar to Appendix D, $(b_{y},b_{z})$ for the correlated Gaussian scenario, denoted by $(b_{y}^{\mathrm{CG}},b_{z}^{\mathrm{CG}})$, is obtained as

[TABLE]

Since $b_{y}^{\mathrm{CG}}$ and $b_{z}^{\mathrm{CG}}$ are again sums of Gaussian RVs, they are also Gaussian RVs. Therefore, with the same reasoning as in Appendix D, $u$ follows a Hoyt (Nakagami-q) distribution with parameters $q=\sqrt{\frac{\min\{\lambda_{1},\lambda_{2}\}}{\max\{\lambda_{1},\lambda_{2}\}}}$ and $\Omega=\lambda_{1}+\lambda_{2}$. Here, $\lambda_{1}$ and $\lambda_{2}$ are now the eigenvalues of $\boldsymbol{\Sigma}_{\mathrm{T}}$ given in Theorem 3. Therefore, formally the same expression for the PDF of the GML is obtained as for independent Gaussian fluctuations. This completes the proof.

Appendix F

In this case, we obtain $\boldsymbol{\Sigma}_{\mathrm{T}}=\boldsymbol{\Sigma}_{\mathrm{CG}}$ which has one non-zero eigenvalue $\lambda_{1}=\zeta^{2}(c_{6}^{2}+c_{7}^{2})$, i.e., $q=0$, and $(b_{y}^{\mathrm{CG}},b_{z}^{\mathrm{CG}})=(c_{6}\delta^{\mathrm{G}},c_{7}\delta^{\mathrm{G}})$ which results in $u=\sqrt{c_{6}^{2}+c_{7}^{2}}\,\,|\delta^{\mathrm{G}}|$, cf. (E) in Appendix E. Since $\delta^{\mathrm{G}}$ follows the Gaussian distribution, $u$ follows a single-sided Gaussian distribution with the PDF given by $f_{u}(u)=\frac{\sqrt{2}}{\sqrt{\pi\lambda_{1}}}\exp\left(-\frac{u^{2}}{2\lambda_{1}}\right)$. Based on the distribution of $u$, the PDF of $f_{h_{g}}(h)$ in (29) is obtained using (21), which completes the proof.

Appendix G

We replace $\epsilon_{s}$ in (D) with $\epsilon_{s}^{\mathrm{CU}},\,\,s\in\{x,y,z,\theta,\phi\}$, cf. (30). Therefore, we have

[TABLE]

where constants $c_{6}$ and $c_{7}$ are given in Theorem 3. Using (55), $u$ is obtained as

[TABLE]

Substituting the uniform distribution in (56) into (21) leads to (31) in Theorem 4. This completes the proof.

Appendix H

The symmetric difference of the first-order Marcum Q-function is approximated by $\underset{(a,b)\to\infty}{\lim}Q(a,b)-Q(b,a)=1-\left(\sqrt{\frac{a}{b}}+\sqrt{\frac{b}{a}}\right)Q\left(a-b\right)$ [37]. Hence, $P_{\mathrm{out}}$ in (33) is simplified to $P_{\mathrm{out}}=\left(\sqrt{\frac{a}{b}}+\sqrt{\frac{b}{a}}\right)$ $Q\left(a-b\right)$. Next, we use $Q(x)\approx\frac{e^{-\frac{1}{2}x^{2}}}{\sqrt{2\pi}x}$ to approximate the Gaussian Q-function at large values, which leads to the simplified expression for the outage probability in (36). This completes the proof.

Acknowledgements

The authors would like to thank Prof. Mohamed-Slim Alouini and Prof. Steve Hranilovic for helpful and insightful discussions on the uniform distribution of the fluctuations of the position and orientation of the UAV during IEEE ICC 2018.