A Geometry-Based Beamforming Channel Model for UAV mmWave Communications

Abstract

Considering the three-dimensional (3D) trajectory, 3D antenna array, and 3D beamforming of unmanned aerial vehicle (UAV), a novel non-stationary millimeter wave (mmWave) geometry-based stochastic model for UAV to vehicle communication channels is proposed. Based on the analysis results of measured and ray tracing simulation data of UAV mmWave communication links, the proposed parametric channel model is constructed by a line-of-sight path, a ground specular path, and two strongest single-bounce paths. Meanwhile, a new parameter computation method is also developed, which is divided into the deterministic (or geometry-based) part and the random (or empirical) part. The simulated power delay profile and power angle profile demonstrate that the statistical properties of proposed channel model are time-variant with respect to the scattering scenarios, positions and beam direction. Moreover, the simulation results of autocorrelation functions fit well with the theoretical ones as well as the measured ones.

1. Introduction

Unmanned aerial vehicles (UAVs) have been expected to be a promising platform to expand the communication range as flying relays in the fifth-generation (5G) and beyond fifth-generation (B5G) communication systems [1,2]. Meanwhile, the millimeter wave (mmWave) from 6 GHz to 100 GHz is also considered to be an important technique to provide high transmission rate and can be applied on UAVs for the small size of mmWave antennas [3,4]. However, the communication distance of UAVs is greatly limited due to the high path loss of mmWave propagation, and thus the three-dimensional (3D) beam-forming technology with antenna arrays has gathered more and more attention [4,5]. Beside these, UAV communication scenarios have some other special characteristics, i.e., 3D flying trajectory and 3D scattering space, which results in the existing channel models for traditional communication scenarios being not suitable any more [6,7]. Therefore, it is essential to deeply understand these new features and obtain a tailored channel model for the UAV mmWave communication system.

There has been growing interest in the UAV channel modeling, which can be classified into deterministic models [8,9] and stochastic models [10,11,12]. Among them, the geometry-based stochastic models (GBSMs) have been widely applied due to the good tradeoff between the generality, accuracy, and complexity [13,14,15,16]. In [13], a 3D non-stationary GBSM for the air-to-ground communication was proposed and both the air and ground terminals were moving. The authors in [14] proposed a novel UAV channel model with two cylinder scatterers around the transceiver. The authors in [15,16] upgraded the 3D GBSM by taking into account of the 3D flying trajectory and arbitrary velocity of UAV. However, all of above GBSMs only focused on the sub-mmWave frequency band and were inapplicable for the UAV mmWave channel.

So far, most of mmWave channel models were mainly aimed at the land mobile communication scenarios [17,18], e.g., tunnel [19,20], microcellular environment [21], and high-voltage substation [22], etc. A few mmWave channel models involving UAV scenarios can be addressed in [23,24,25,26]. In [23,24], the authors used ray tracing (RT) simulated data to develop UAV mmWave channel models and analyze the channel parameters, i.e., received power, path delay, angle, etc. It should be noted that the channel model was 2D in [23] and the ground terminal was fixed in [24]. Moreover, the RT-based channel modeling is generally time-consuming and has high complexity. In [25], a 3D non-stationary mmWave channel for UAV communication was proposed based on the geometry-based stochastic model (GBSM) method, but the flight velocity was constant and the receiving terminal was also fixed. Recently, the authors in [26] proposed a mmWave UAV channel model allowing 3D trajectories, but the rotation of 3D-shaped antenna array and the effect of beam-forming were not considered.

The GBSM is a great alternative to model the UAV mmWave beam channel, but we note that the existing GBSMs for UAV mmWave channel cannot completely cover the new features of UAV to vehicle (U2V) beam communications. Especially, the characteristic of 3D beam-forming is not included in the related literature, which would change the covering range of signal or the distribution of scatterers and furthermore affects the channel properties, i.e., path angles and power gains. This paper aims to fill these research gaps. The major contributions and novelties of this paper are summarized as follows:

(1)

A 3D GBSM for U2V mmWave beam channel considering the 3D arbitrary trajectory, 3D antenna array, and 3D beam-forming of UAV is proposed. To achieve the tradeoff between generality, accuracy, and complexity, the model only takes into account the line-of-sight (LoS) path, ground specular (GS) path, and two strongest single-bounce (SB) paths.

(2)

A hybrid computation method of channel parameters, i.e., geometry-based parameters and data-based parameters, for the proposed model is developed. The geometry-based parameters, e.g., the locations of terminals, the mean angles and delays of paths, are calculated by the time-variant but deterministic geometric relationships, and the data-based channel parameters, e.g., the angle offset and delay offset of the rays, the path powers, are generated randomly from the corresponding distribution fitted by RT simulation or measured data.

(3)

Considering an urban U2V mmWave communication scenario, the channel parameters, i.e., path delays, received powers, and angles, are simulated and demonstrated. Moreover, the simulation results of the second order statistical properties, i.e., autocorrelation function (ACF) and Doppler power spectral density (DPSD), are also validated by theoretical and measured ones.

The rest paper is organized as follows. In Section 2, a 3D GBSM for U2V mmWave beam channel is proposed. Section 3 gives the computation method of geometry-based parameters and data-based parameters. The simulation and analytical results of the channel parameters and statistical properties are given in Section 4. Finally, conclusions are drawn in Section 5.

2. UAV mmWave Channel Model

Let us consider a typical U2V communication scenario as shown in Figure 1, where the 3D beamforming is applied to compensate the high path loss. The vehicle is equipped with omnidirectional antennas. Within the beam, the gain coefficient from the antenna array varies with the angle offset between the beam center and different paths, which would affect the power gain of propagation channel. In Figure 1, two independent coordinate systems are denoted as the UAV coordinate system and vehicle coordinate system with their origins at the central of UAV and vehicle, respectively, and the arbitrary velocities of UAV or vehicle can be denoted as

$${\mathbf{v}}_{\mathrm{T}/\mathrm{R}}\left(t\right)=v_{\mathrm{T}/\mathrm{R}}\left(t\right)\left[\begin{array}{c}cos{\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)cos{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)\\ sin{\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)cos{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)\\ sin{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)\end{array}\right]$$

(1)

where $v_{\mathrm{T}/\mathrm{R}}\left(t\right)$, ${\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)$ and ${\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)$ are the magnitude, azimuth angle, and elevation angle of ${\mathbf{v}}_{\mathrm{T}/\mathrm{R}}\left(t\right)$, respectively. It should be noted that ${\left(·\right)}_{\mathrm{T}/\mathrm{R}}$ represents two equations for T and R, respectively.

The U2V beamforming propagation channel normally includes a LoS path and tremendous NLoS paths, i.e., GS path, SB paths, double-bounce (DB) paths, etc. To simplify the channel model, we have performed large amount of RT simulations and obtained massive raw data of mmWave U2V channels in 28 GHz under four scenarios, i.e., urban, forest, hill, and sea. It is assumed that the LoS path always exists during the simulation. As shown in Figure 2a, we find that the powers of LoS path and three strongest paths are at least $90\%$ and $95\%$ of the total one, respectively. Moreover, the summed power of four strongest paths is over $99\%$ of the total one and thus the path number of channel model can be simplified into 4. Furthermore, Figure 2b gives more details of the received powers of different paths. In the figure, the UAV flies through several different trajectories under the urban scenario and the received powers are averaged. It can be seen that the power of the LoS path is normally 20 dB and 40 dB more than the GS path and the SB paths, respectively. Moreover, the power of the strongest DB path is lower than the ones of strongest SB path and second strongest SB path, and thus the four strongest paths normally represent the LoS path, GS path, 1st SB path, and 2st SB path.

Moreover, it is assumed that the LoS path includes one direct ray and each NLoS path includes several non-direct rays with similar delays. Therefore, in this paper the channel impulse response (CIR) between the pth UAV antenna and the qth vehicle antenna scaled by the K-factor K only includes the LoS path $h_{}^{\mathrm{LoS}}\left(t\right)$ and three strongest NLoS paths, i.e., $h_{}^{\mathrm{GS}}(\tau ,t)$, $h_{}^{\mathrm{S}{\mathrm{B}}_1}(\tau ,t)$, $h_{}^{\mathrm{S}{\mathrm{B}}_2}(\tau ,t)$ as

$$\begin{array}{cc}\hfill h(\tau ,t)& =h_{}^{\mathrm{LoS}}\left(t\right)+\underset{groundspecular}{\underbrace{h_{}^{\mathrm{GS}}(\tau ,t)}}+\underset{Singlebounce}{\underbrace{h_{}^{\mathrm{S}{\mathrm{B}}_1}(\tau ,t)+h_{}^{\mathrm{S}{\mathrm{B}}_2}(\tau ,t)}}\hfill \\ \hfill & =\sqrt{\frac{K\left(t\right)}{K\left(t\right)+1}}A({\theta }^{\mathrm{LoS}},{\phi }^{\mathrm{LoS}}){\tilde{h}}_{}^{\mathrm{LoS}}\left(t\right)\hfill \\ \hfill & +\sqrt{\frac{1}{K\left(t\right)+1}}\sum _{j=1}^3\sum _{m=1}^{M}A({\theta }_m^j,{\phi }_m^j)\sqrt{P_m^j\left(t\right)}{\tilde{h}}_m^j\left(t\right)\delta (t-{\tau }_m^j)\hfill \end{array}$$

(2)

where M is the ray number of each NLoS path, $j\in \left\{1,2,3\right\}$ represents the GS path, SB1 path and SB2 path, respectively. In (2), $A({\theta }^{\mathrm{LoS}},{\phi }^{\mathrm{LoS}})$ and $A({\theta }_m^j,{\phi }_m^j)$ are the gain coefficient of LoS path and NLoS paths within the beam, respectively, $P_m^j\left(t\right)$ and ${\tau }_m^j\left(t\right)$ are the power and delay of the mth ray within the jth NLoS path, respectively. Moreover, ${\tilde{h}}_m^j\left(t\right)$ denotes the channel coefficient of mth ray and can be modeled as

$$\begin{array}{cc}\hfill {\tilde{h}}_m^j\left(t\right)& =exp\left(j{\mathsf{\Phi }}_{\mathrm{I}}\right)exp\left(\mathrm{j}2\pi \frac{{\mathbf{r}}_{\mathrm{R},m}^j\left(t\right)·{\mathbf{R}}_{\mathrm{R}}^{}\left(t\right)·{\mathbf{d}}_{\mathrm{R}}\left(t\right)}{\lambda }\right)\hfill \\ & ·exp\left(\mathrm{j}2\pi \frac{{\mathbf{r}}_{\mathrm{T},m}^j\left(t\right)·{\mathbf{R}}_{\mathrm{T}}\left(t\right)·{\mathbf{d}}_{\mathrm{T}}\left(t\right)}{\lambda }\right)·exp\left(\mathrm{j}2\pi {\int }_{t_0}^t{f}_m^j\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }\right)\hfill \end{array}$$

(3)

where ${\mathsf{\Phi }}_{\mathrm{I}}$ represents the random initial phase distributed uniformly over $[0,2\pi )$, $\lambda$ is the wavelength, $t_0$ is the initial time instant, and ${\mathbf{d}}_{\mathrm{T}/\mathrm{R}}\left(t\right)$ denotes the location vectors of UAV transmitting antenna (or vehicle receiving antenna) in their own coordinate systems and can be described as

$${\mathbf{d}}_{\mathrm{T}/\mathrm{R}}\left(t\right)=\left[\begin{array}{c}{d}_{\mathrm{T}/\mathrm{R}}^x\left(t\right)\\ d_{\mathrm{T}/\mathrm{R}}^y\left(t\right)\\ d_{\mathrm{T}/\mathrm{R}}^z\left(t\right)\end{array}\right]$$

(4)

where $d_{\mathrm{T}/\mathrm{R}}^x\left(t\right)$, $d_{\mathrm{T}/\mathrm{R}}^y\left(t\right)$, and $d_{\mathrm{T}/\mathrm{R}}^z\left(t\right)$ represent the x, y, and z component of ${\mathbf{d}}_{\mathrm{T}/\mathrm{R}}\left(t\right)$. In (3), ${\mathbf{r}}_{\mathrm{T},m}^j\left(t\right)$ (or ${\mathbf{r}}_{\mathrm{R},m}^j\left(t\right)$) and $f_m^j$ are the spherical unit vectors and Doppler frequency of the mth ray, respectively, and ${\mathbf{r}}_{\mathrm{T},m}^j\left(t\right)$ (or ${\mathbf{r}}_{\mathrm{R},m}^j\left(t\right)$) can be further expressed as

$${\mathbf{r}}_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)=\left[\begin{array}{c}cos{\beta }_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)cos{\alpha }_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)\\ cos{\beta }_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)sin{\alpha }_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)\\ sin{\beta }_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)\end{array}\right]$$

(5)

where ${\alpha }_{\mathrm{T}/\mathrm{R},m}^j$ is the azimuth angle of departure (AAoD) or arrival (AAoA), ${\beta }_{\mathrm{T}/\mathrm{R},m}^j$ is the elevation angle of departure (EAoD) or arrival (EAoA). During the movement of UAV and vehicle, the location of each antenna may change, and in this paper a rotation matrix ${\mathbf{R}}_{\mathrm{T}/\mathrm{R}}\left(t\right)$ is introduced to take this factor into account as

$${\mathbf{R}}_{\mathrm{T}/\mathrm{R}}\left(t\right)=\left[\begin{array}{ccc}cos{\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)cos{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)& -sin{\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)& -cos{\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)sin{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)\\ sin{\alpha }_{\mathrm{tx}/\mathrm{rx}}^v\left(t\right)cos{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)& cos{\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)& -sin{\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)sin{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)\\ sin{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)& 0& cos{\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)\end{array}\right]$$

(6)

The LoS path between the UAV and vehicle can be viewed as a special case of NLoS path and the corresponding channel coefficient can be expressed as

$$\begin{array}{c}{\tilde{h}}_{}^{\mathrm{LoS}}\left(t\right)=exp\left(-\mathrm{j}2\pi \frac{D_{\mathrm{LoS}}\left(t\right)}{\lambda }\right)exp\left(\mathrm{j}2\pi \frac{{\mathbf{r}}_{\mathrm{R}}^{\mathrm{LoS}}\left(t\right)·{\mathbf{R}}_{\mathrm{R}}\left(t\right)·{\mathbf{d}}_{\mathrm{R}}\left(t\right)}{\lambda }\right)\hfill \\ ·exp\left(\mathrm{j}2\pi \frac{{\mathbf{r}}_{\mathrm{T}}^{\mathrm{LoS}}\left(t\right)·{\mathbf{R}}_{\mathrm{T}}\left(t\right)·{\mathbf{d}}_{\mathrm{T}}\left(t\right)}{\lambda }\right)exp\left(\frac{\mathrm{j}2\pi }{\lambda }{\int }_0^t{f}_{}^{\mathrm{LoS}}\left(t^{\prime }\right))d{t}^{\prime }\right)\delta (t-{\tau }^{\mathrm{LoS}}\left(t\right))\hfill \end{array}$$

(7)

where $D_{\mathrm{LoS}}\left(t\right)$ is the distance between the UAV and vehicle, ${\mathbf{r}}_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)$ and $f_{}^{\mathrm{LoS}}$ denote the spherical unit vectors and Doppler frequency of LoS path, respectively, and ${\mathbf{r}}_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)$ can be expressed by ${\alpha }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}$ and ${\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}$ as

$$\begin{array}{c}{\mathbf{r}}_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)=\left[\begin{array}{c}cos{\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)cos{\alpha }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)\\ cos{\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)sin{\alpha }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)\\ sin{\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)\end{array}\right].\hfill \end{array}$$

(8)

3. Hybrid Computation Method of Channel Parameters

3.1. Geometry-Based Parameters

The deterministic part of channel parameters can be calculated according to the geometric relationships of communication scenario, i.e., locations and velocities of transceivers and scatterers. Since the UAV and vehicle move along with 3D arbitrary trajectories, the time-variant location vector of UAV (or vehicle) can be expressed as

$${\mathbf{L}}_{\mathrm{T}}\left(t\right)={\mathbf{L}}_{\mathrm{T}}\left(t_0\right)+{\int }_{t_0}^t{\mathbf{v}}_{\mathrm{T}}\left(t\right)\mathrm{d}t$$

(9)

$${\mathbf{L}}_{\mathrm{R}}\left(t\right)={\mathbf{L}}_{\mathrm{R}}\left(t_0\right)+{\int }_{t_0}^t{\mathbf{v}}_{\mathrm{R}}\left(t\right)\mathrm{d}t$$

(10)

where ${\mathbf{L}}_{\mathrm{T}/\mathrm{R}}\left(t_0\right)$ denotes the initial location vector of UAV (or vehicle) at $t=t_0$. The distance vector between the UAV and vehicle ${\mathbf{D}}_{\mathrm{LoS}}^{}\left(t\right)$ can be expressed as

$$\begin{array}{cc}\hfill {\mathbf{D}}_{\mathrm{LoS}}^{}\left(t\right)& ={\mathbf{L}}_{\mathrm{T}}^{}\left(t\right)-{\mathbf{L}}_{\mathrm{R}}^{}\left(t\right)\hfill \\ \hfill & =D_{\mathrm{LoS}}^{}\left(t_0\right){\mathbf{r}}_{\mathrm{T}}^{\mathrm{LoS}}\left(t\right)+{\int }_{t_0}^t{\mathbf{v}}_{\mathrm{T},\mathrm{R}}\left(t\right)\mathrm{d}t\hfill \end{array}$$

(11)

where $D_{\mathrm{LoS}}^{}\left(t_0\right)$ is the initial distance of LoS path and can be expressed as

$$D_{\mathrm{LoS}}^{}\left(t_0\right)=∥{\mathbf{L}}_{\mathrm{T}}\left(t_0\right)-{\mathbf{L}}_{\mathrm{R}}\left(t_0\right)∥$$

(12)

And ${\mathbf{v}}_{\mathrm{T},\mathrm{R}}\left(t\right)$ is the relative velocity between the UAV and vehicle which can be expressed as

$${\mathbf{v}}_{\mathrm{T},\mathrm{R}}\left(t\right)={\mathbf{v}}_{\mathrm{T}}\left(t\right)-{\mathbf{v}}_{\mathrm{R}}\left(t\right).$$

(13)

Similarly, the distance vector between the UAV (or vehicle) and jth scatterer ${\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^{}\left(t\right)$ can be expressed as

$$\begin{array}{cc}\hfill {\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^{}\left(t\right)& ={\mathbf{L}}_{\mathrm{T}/\mathrm{R}}^{}\left(t\right)-{\mathbf{L}}_j^{}\left(t\right)\hfill \\ \hfill & =D_{\mathrm{T}/\mathrm{R},j}^{}\left(t_0;t\right){\mathbf{r}}_{\mathrm{T}/\mathrm{R}}^j\left(t\right)+{\int }_{t_0}^t{\mathbf{v}}_{\mathrm{T}/\mathrm{R}}\left(t\right)\mathrm{d}t\hfill \end{array}$$

(14)

where $D_{\mathrm{T}/\mathrm{R},j}^{}\left(t_0;t\right)$ is the initial distance between the UAV (or vehicle) and $j\mathrm{th}$ scatterer. In (14), ${\mathbf{r}}_{\mathrm{T}/\mathrm{R}}^j\left(t\right)$ is the spherical unit vectors of each NLoS path, which can be expressed by the mean angles denoted by ${\overline{\alpha }}_{\mathrm{T}/\mathrm{R}}^j$ and ${\overline{\beta }}_{\mathrm{T}/\mathrm{R}}^j$ as

$${\mathbf{r}}_{\mathrm{T}/\mathrm{R}}^j\left(t\right)=\left[\begin{array}{c}cos{\overline{\beta }}_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)cos{\overline{\alpha }}_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)\\ cos{\overline{\beta }}_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)sin{\overline{\alpha }}_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)\\ sin{\overline{\beta }}_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)\end{array}\right]$$

(15)

Moreover, the scatterers in this paper are assumed to be static and the valid ones covered by the 3D beam will change along with the beam direction. To evolve the geometry parameter of distance, the initial distance should be refreshed when the distribution of the valid scatterers have changed. Therefore, the initial distance $D_{\mathrm{T}/\mathrm{R},j}^{}\left(t_0;t\right)$ can be rewritten as

$$D_{\mathrm{T}/\mathrm{R},j}^{}\left(t_0;t\right)=D_{\mathrm{T}/\mathrm{R},j}^{}(t_0+k{T}_0)\mathrm{W}\left(t-k{T}_0\right),k=0,1,2...$$

(16)

where the window function $\mathrm{W}(·)$ is introduced and can be denoted as

$$\mathrm{W}\left(t\right)\stackrel{\mathsf{\Delta }}{=}\left\{\begin{array}{c}1{,}_{}^{}0\le t\le T_0\\ 0{,}_{}^{}\mathrm{otherwise}\end{array}\right.$$

(17)

where $T_0$ is the stationary interval of scatterers and it is related with the beam width and the velocity of both terminals. The distance between UAV and vehicle in the LoS scenario and the one between the UAV (or vehicle) and jth scatterer can be calculated respectively by

$$D_{\mathrm{LoS}}^{}\left(t\right)=\sqrt{\begin{array}{c}{\left(D_{\mathrm{LoS}}^{}\left(t_0\right)cos({\alpha }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t_0\right))cos({\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t_0\right))+{\int }_{t_0}^t\left(v_{\mathrm{T},\mathrm{R}}\left(t\right)\right)·cos({\alpha }_{\mathrm{T},\mathrm{R}}^v\left(t\right))·cos({\beta }_{\mathrm{T},\mathrm{R}}^v\left(t\right))dt\right)}^2\hfill \\ +{\left(D_{\mathrm{LoS}}^{}\left(t_0\right)cos({\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t_0\right))sin({\alpha }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t_0\right))+{\int }_{t_0}^t\left(v_{\mathrm{T},\mathrm{R}}\left(t\right)\right)·sin({\alpha }_{\mathrm{T},\mathrm{R}}^v\left(t\right))cos({\beta }_{\mathrm{T},\mathrm{R}}^v\left(t\right))dt\right)}^2\hfill \\ +{\left(D_{\mathrm{LoS}}^{}\left(t_0\right)sin({\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t_0\right))+{\int }_{t_0}^t\left(v_{\mathrm{T},\mathrm{R}}\left(t\right)\right)·sin({\beta }_{\mathrm{T},\mathrm{R}}^v\left(t\right))dt\right)}^2\hfill \end{array}}$$

(18)

$$D_{\mathrm{T}/\mathrm{R},j}^{}\left(t\right)=\sqrt{\begin{array}{c}{\left(D_{\mathrm{T}/\mathrm{R},j}^{}\left(t_0;t\right)cos({\overline{\alpha }}_{\mathrm{T}/\mathrm{R}}^j\left(t_0\right))cos({\overline{\beta }}_{\mathrm{T}/\mathrm{R}}^j\left(t_0\right))+{\int }_{t_0}^t{v}_{\mathrm{T}/\mathrm{R}}\left(t\right)cos({\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right))cos({\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right))dt\right)}^2\hfill \\ +{\left(D_{\mathrm{T}/\mathrm{R},j}^{}\left(t_0;t\right)cos({\overline{\beta }}_{\mathrm{T}/\mathrm{R}}^j\left(t_0\right))sin({\overline{\alpha }}_{\mathrm{T}/\mathrm{R}}^j\left(t_0\right))+{\int }_{t_0}^t{v}_{\mathrm{T}/\mathrm{R}}\left(t\right)sin({\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right))cos({\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right))dt\right)}^2\hfill \\ +{\left(D_{\mathrm{T}/\mathrm{R},j}^{}\left(t_0;t\right)sin({\overline{\beta }}_{\mathrm{T}/\mathrm{R}}^j\left(t_0\right))+{\int }_{t_0}^t{v}_{\mathrm{T}/\mathrm{R}}\left(t\right)sin({\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right))dt\right)}^2\hfill \end{array}}$$

(19)

where ${\alpha }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)$ and ${\beta }_{\mathrm{T}/\mathrm{R}}^v\left(t\right)$ denote the relative travel direction between the UAV and vehicle on the azimuth and elevation plane, respectively.

Based on the geometric relationships, the time-variant angles such as the EAoD, AAoD, EAoA, and AAoA of LoS path under dynamic U2V communication scenarios can be obtained respectively by

$${\alpha }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)=\left\{\begin{array}{c}arccos(\frac{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^x\left(t\right)∥}{\sqrt{{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^x\left(t\right)∥}^2+{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^y\left(t\right)∥}^2}}),{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^x\left(t\right)\ge 0\hfill \\ \pi -arccos(\frac{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^x\left(t\right)∥}{\sqrt{{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^x\left(t\right)∥}^2+{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^y\left(t\right)∥}^2}}),{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^x\left(t\right)<0\hfill \end{array}\right.$$

(20)

$${\beta }_{\mathrm{T}/\mathrm{R}}^{\mathrm{LoS}}\left(t\right)=arcsin(\frac{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^z\left(t\right)∥}{D_{\mathrm{LoS}}^{}\left(t\right)})$$

(21)

where ${\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^x\left(t\right)$, ${\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^y\left(t\right)$, ${\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^z\left(t\right)$ denote the x, y, z component of ${\mathbf{D}}_{\mathrm{T}/\mathrm{R},\mathrm{LoS}}^{}\left(t\right)$, respectively. Then, the Doppler frequency of LoS path can be rewritten as

$$\begin{array}{cc}\hfill f_{}^{\mathrm{LoS}}\left(t\right)& =\frac{{\mathbf{v}}_{\mathrm{T}}\left(t\right)·{\mathbf{r}}_{\mathrm{T}}^{\mathrm{LoS}}\left(t\right)}{\lambda }+\frac{{\mathbf{v}}_{\mathrm{R}}\left(t\right)·{\mathbf{r}}_{\mathrm{R}}^{\mathrm{LoS}}\left(t\right)}{\lambda }\hfill \\ \hfill & =\frac{v_{\mathrm{T}}\left(t\right)\left(cos({\alpha }_{\mathrm{T}}^{\mathrm{LoS}}\left(t\right)-{\alpha }_{\mathrm{T}}^v\left(t\right))cos{\beta }_{\mathrm{T}}^{\mathrm{LoS}}\left(t\right)cos{\beta }_{\mathrm{T}}^v\left(t\right)+sin{\beta }_{\mathrm{T}}^{\mathrm{LoS}}\left(t\right)sin{\beta }_{\mathrm{T}}^v\left(t\right)\right)}{\lambda }\hfill \\ \hfill & +\frac{+v_{\mathrm{R}}\left(t\right)\left(cos({\alpha }_{\mathrm{R}}^{\mathrm{LoS}}\left(t\right)-{\alpha }_{\mathrm{R}}^v\left(t\right))\cos {\beta }_{\mathrm{R}}^{\mathrm{LoS}}\left(t\right)cos{\beta }_{\mathrm{R}}^v\left(t\right)+sin{\beta }_{\mathrm{R}}^{\mathrm{LoS}}\left(t\right)sin{\beta }_{\mathrm{R}}^v\left(t\right)\right)}{\lambda }.\hfill \end{array}$$

(22)

For the NLoS paths, the mean angles of time-variant EAoD, AAoD or EAoA, AAoA can be calculated respectively by

$${\overline{\alpha }}_{\mathrm{T}/\mathrm{R}}^j\left(t\right)=\left\{\begin{array}{c}arccos(\frac{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^x\left(t\right)∥}{\sqrt{{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^x\left(t\right)∥}^2+{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^y\left(t\right)∥}^2}}),{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^x\left(t\right)\ge 0\hfill \\ \pi -arccos(\frac{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^x\left(t\right)∥}{\sqrt{{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^x\left(t\right)∥}^2+{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^y\left(t\right)∥}^2}}),{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^x\left(t\right)<0\hfill \end{array}\right.$$

(23)

$${\overline{\beta }}_{\mathrm{T}/\mathrm{R}}^j\left(t\right)=arcsin(\frac{∥{\mathbf{D}}_{\mathrm{T}/\mathrm{R},j}^z\left(t\right)∥}{D_{\mathrm{T}/\mathrm{R},j}^{}\left(t\right)}).$$

(24)

Similarly, the time-variant delays of LoS and NLoS paths are related with the transmission distance and they can be calculated respectively by

$${\tau }_{}^{\mathrm{LoS}}\left(t\right)=\frac{D_{\mathrm{LoS}}^{}\left(t\right)}{c}$$

(25)

$${\overline{\tau }}_{}^j\left(t\right)=\frac{D_{\mathrm{T},j}^{}\left(t\right)+D_{\mathrm{R},j}^{}\left(t\right)}{c}$$

(26)

where c is the speed of light.

3.2. Data-Based Stochastic Parameters

Note that the channel parameters obtained by the geometric relationships cannot reflect the statistical properties of random rays within the NLoS path. In this paper, we have conducted simulations in 28 GHz based on the RT method for up to 7200 different UAV channels under four scenarios, i.e., urban, forest, hill, and sea. Huge amount of obtained channel data is used to analyze the stochastic characteristic of random rays. Based on these analytical results, the computation method of intra-path or ray parameters are developed as follows.

Each NLoS path may include several rays with different delays, which can be described as the relative delay between the intra-path ray delay and mean path delay. The data-based analytical results of ray delay offset under different scenarios are shown in Figure 3. As we can see that it is more likely to arise the delay offset under the urban scenario and the value of the ray delay offset with high probability can be up to 10 ns since the urban scenario has rich scatterers. In the other simulation scenarios, the values of the ray delay offset with high probability are normally less than 6 ns.

Moreover, it can be seen from Figure 3 that the PDFs of delay offset fit well with modified Gaussian distribution with zero mean value. Therefore, the delay offset of the mth ray within jth path can be modeled as

$$f(\mathsf{\Delta }{\tau }_m^j)=\frac{a_{\tau }}{\sqrt{2\pi }{\sigma }_{\tau }}exp(-\frac{{(\mathsf{\Delta }{\tau }_m^j)}^2}{2{{\sigma }_{\tau }}^2})+b_{\tau }$$

(27)

where $a_{\tau }$ and $b_{\tau }$ are the weight factor and correction factor of the magnitude respectively, and ${\sigma }_{\tau }^{}$ is the standard deviation of delay offset. In Figure 3, $a_{\tau }$ is $9.36$, $9.80$, $9.56$, $10.20$, ${\sigma }_{\tau }^{}$ is $11.24$, $5.10$, $8.11$, $8.13$ and $b_{\tau }$ is $0.011$, $0.005$, $0.015$, $0.016$ in four different scenarios. Then, the ray delay can be generated by

$${\tau }_m^j\left(t\right)={\overline{\tau }}_{}^j\left(t\right)+\mathsf{\Delta }{\tau }_m^j.$$

(28)

Furthermore, the power of each ray can be calculated according to the ray delay by the exponential distribution as

$$P_m^j\left(t\right)=\xi exp\left(-1000\xi {\tau }_m^j\left(t\right)\right)$$

(29)

where $\xi$ is the rate parameter of exponential distribution. According to the analytical results of the simulation data, $\xi$ is well fitted by $5.85$, $26.7$, $22.8$, and $25.05$ in four typical scenarios. When the power of LoS path is normalized to be 0 dB, the power of each ray can be normalized as

$${\tilde{P}}_m^j\left(t\right)=\frac{P_m^j\left(t\right)}{A({\theta }^{\mathrm{LoS}},{\phi }^{\mathrm{LoS}})·K\left(t\right)·{\displaystyle \sum _{m=1}^M}{P}_m^j\left(t\right)}$$

(30)

where the total cluster powers should be $1/\left(K\right(t)+1)$ as shown in (2).

Similarly, the mean angles of the NLoS path are along with ray angle offset as well. The analytical results of the azimuth and elevation angle offset are shown in Figure 4. As can be seen that the values of the azimuth angle offset with high probability are smaller than the ones of elevation angle offset, because the width of the scatterers is normally shorter than the height of the scatterers especially under the urban scenario. Moreover, we can find that the PDFs of the azimuth and elevation angle offset can be fitted well by modified Gaussian distribution and Laplacian distribution respectively, which is similar with the recommendation in [27]. For the azimuth offset of AAoA or AAoD, the PDFs can be described by $a_{\alpha }$, ${\sigma }_{\alpha }^{}$, and $\mathsf{\Delta }{\alpha }_m^j$ as

$$f(\mathsf{\Delta }{\alpha }_m^j)=\frac{a_{\alpha }}{\sqrt{2\pi }{\sigma }_{\alpha }}exp(-\frac{{(\mathsf{\Delta }{\alpha }_m^j)}^2}{2{{\sigma }_{\alpha }}^2})+b_{\alpha }.$$

(31)

In Figure 4a, $a_{\alpha }$ is $0.51$, $0.98$, $0.68$, $0.89$, ${\sigma }_{\alpha }^{}$ is $0.81$, $0.51$, $0.60$, $0.60$ and $b_{\alpha }$ is $0.017$, $0.006$, $0.014$, $0.013$, respectively. The PDFs of elevation offset EAoA or EAoD can be expressed by $\mathsf{\Delta }{\beta }_m^j$ as

$$f(\mathsf{\Delta }{\beta }_m^j)=\frac{a_{\beta }}{2{\sigma }_{\beta }}exp(-\frac{\left|\mathsf{\Delta }{\beta }_m^j\right|}{{\sigma }_{\beta }}).$$

(32)

In Figure 4b, $a_{\beta }$ is $1.810$, $1.008$, $1.005$, $1.007$ and ${\sigma }_{\beta }$ is $3.52$, $0.57$, $0.70$, $0.58$ in different scenarios. Then the ray angle can be generated by

$${\alpha }_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)={\overline{\alpha }}_{\mathrm{T}/\mathrm{R}}^j\left(t\right)+\mathsf{\Delta }{\alpha }_m^j$$

(33)

$${\beta }_{\mathrm{T}/\mathrm{R},m}^j\left(t\right)={\overline{\beta }}_{\mathrm{T}/\mathrm{R}}^j\left(t\right)+\mathsf{\Delta }{\beta }_m^j.$$

(34)

Finally, the Doppler frequency of NLoS path should be calculated by the relative velocity both between the UAV and scatterers and between the vehicle and scatterers. However, the scatterers in this paper are assumed to be static and then the Doppler frequency of NLoS path can be got equivalently by the velocity of the UAV and vehicle. Thus, the Doppler frequency of each ray $f_m^j\left(t\right)$ within jth NLoS path can be described as

$$f_m^j\left(t\right)=f_{\mathrm{T},m}^j\left(t\right)+f_{\mathrm{R},m}^j\left(t\right)$$

(35)

where $f_{\mathrm{T},m}^j\left(t\right)$ and $f_{\mathrm{R},m}^j\left(t\right)$ can be further expressed as

$$\begin{array}{cc}\hfill f_{\mathrm{T},m}^j\left(t\right)& =\frac{{\mathbf{v}}_{\mathrm{T}}\left(t\right)·{\mathbf{r}}_{\mathrm{T},m}^j\left(t\right)}{\lambda }\hfill \\ \hfill & =\frac{v_{\mathrm{T}}\left(t\right)\left(cos({\alpha }_{\mathrm{T},m}^j\left(t\right)-{\alpha }_{\mathrm{T}}^v\left(t\right))\cos {\beta }_{\mathrm{T},m}^j\left(t\right)cos{\beta }_{\mathrm{T}}^v\left(t\right)+sin{\beta }_{\mathrm{T},m}^j\left(t\right)sin{\beta }_{\mathrm{T}}^v\left(t\right)\right)}{\lambda }\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill f_{\mathrm{R},m}^j\left(t\right)& =\frac{{\mathbf{v}}_{\mathrm{R}}\left(t\right)·{\mathbf{r}}_{\mathrm{R},m}^j\left(t\right)}{\lambda }\hfill \\ \hfill & =\frac{v_{\mathrm{R}}\left(t\right)\left(cos({\alpha }_{\mathrm{R},m}^j\left(t\right)-{\alpha }_{\mathrm{R}}^v\left(t\right))\cos {\beta }_{\mathrm{R},m}^j\left(t\right)cos{\beta }_{\mathrm{R}}^v\left(t\right)+sin{\beta }_{\mathrm{R},m}^j\left(t\right)sin{\beta }_{\mathrm{R}}^v\left(t\right)\right)}{\lambda }.\hfill \end{array}$$

(37)

4. Simulation Results and Validations

To illustrate and verify the proposed U2V beamforming channel model, we take the urban scenario as an example. The 3D trajectories of both terminals and 3D beam are considered. It should be mentioned that the angular error of beam tracking is ignored and the rest simulation parameters are shown in

Table 1. Simulation parameters.

Definition	Value	Definition	Value
$v_{\mathrm{T}}\left(t\right)$	$10+0.5t$ m/s	$v_{\mathrm{R}}\left(t\right)$	$2+t$ m/s
${\alpha }_{\mathrm{T}}^v\left(t\right)$	$120-2t^{\circ }$	${\beta }_{\mathrm{T}}^v\left(t\right)$	$6+5t^{\circ }$
${\alpha }_{\mathrm{R}}^v\left(t\right)$	$-120+2t^{\circ }$	${\beta }_{\mathrm{R}}^v\left(t\right)$	$0^{\circ }$
$\lambda$	$3/280$ m	K	7 dB
$∥{\mathbf{L}}_{\mathrm{T}}\left(t_0\right)∥$	400 m	$∥{\mathbf{L}}_{\mathrm{R}}\left(t_0\right)∥$	100 m
${\sigma }_{\tau }^{}$	$11.24$	${\sigma }_{\alpha }^{}$	$0.81$
${\sigma }_{\beta }^{}$	$3.52$	M	12

. Moreover, the generation method of 3D beam in [28] is adopted in this paper and the beam width is assumed to be $22.5$ degrees. The power gain coefficient of 3D beam is simulated and given in Figure 5. As we can see that the power gain will decay with the azimuth angle and elevation angle deviating from the beam center.

The delay, angle, and power of each path and intra-path ray become complicated due to the movement of terminals and 3D beamforming. Based on the above computation method, we run the proposed channel model and obtain the time-variant normalized PDPs at three different time instants $t_1$ = 0 s, $t_2$ = 5 s, and $t_3$ = 10 s as shown in Figure 6. As we can see, each NLoS path includes several intra-path rays. Although the path power and intra-path ray power will be affected a little by the gain coefficient of 3D beam, the trend of mean path power in dB is linear attenuation. It should be noted that the exponential expression in (29) can be rewritten to a linear expression with respect to the power in dB.

Similarly, the angle parameters can be simulated according to the geometric parameters of the simulation scenario and above statistical distribution of the angle offset. The time-variant PAPs of LoS path and rays within three NLoS paths at three different time instants $t_1$ = 0 s, $t_2$ = 5 s, and $t_3$ = 10 s are shown in Figure 7. As we can see from the Figure 7a, taking the AAoD and EAoD of LoS path as the central point, the deviation range of NLoS paths at different time instant are all limited within 22.5 degrees which is consistent with the beam width. In addition, the beam points to different direction at different time instants due to the movement of the UAV and vehicle. In Figure 7b, the spread range of AAoA is up to 200 degrees which is much wider than the EAoA since the scatterers may distribute around the vehicle arbitrary but the height of the scatterers is limited. Moreover, the vehicle antenna is usually close to the ground and thus the EAoA would be within $\left[0^{\circ }{}_{}^{}{90}^{\circ }\right]$.

The ACF is an important second order statistical property which describes the fading correlation over time. It is usually used to verify the correctness of theoretical channel model. The theoretical normalized ACF of proposed model by the definition can be expressed as

$$\begin{array}{cc}\hfill r_{hh}\left(\mathsf{\Delta }t;t\right)& =\frac{\mathrm{E}\left[h^{*}\left(t\right)h(t+\mathsf{\Delta }t)\right]}{\sqrt{\mathrm{E}\left[{\left|h^{*}\left(t\right)\right|}^2\right]\mathrm{E}\left[{\left|h(t+\mathsf{\Delta }t)\right|}^2\right]}}\hfill \end{array}$$

(38)

where $\mathrm{E}\left(·\right)$ is the expectation function and ${\left(·\right)}^{*}$ is the complex conjugate.

By submitting $\left(2\right)$ into $\left(38\right)$, we can obtain the theoretical ACF of proposed model. On the other hand, the simulated ACF can be obtained by using the generated channel data. The theoretical and simulated ACFs are compared in Figure 8 at three different time instants $t_1$ = 0 s, $t_2$ = 5 s, and $t_3$ = 10 s. It can be seen that the ACFs change over time due to the time-variant channel parameters and the simulated results fit well with the theoretical ones. Furthermore, we can get the DPSDs by using the Fourier transform of ACFs and the simulated results of the 3D time-variant DPSDs are shown in Figure 9. As was shown, the Doppler frequency changes in a complicated way under the U2V communication scenario due to the 3D trajectory and 3D beamforming. It should be noted that all the LoS path and NLoS paths are considered together in the simulation results of the ACFs and DPSDs.

To further verify the consistency of proposed channel model with realistic channels, the simulated ACF is compared with the measured one in [29]. It should be mentioned that there is very little literature involving UAV mmWave channel measurement [30,31,32] and none of them so far analyzed the measured ACFs. Noting that GBSM is a general channel modeling method which can be applied in both the sub-mmWave band and mmWave band, the proposed model could be used for the sub-mmWave GBSM by adjusting some channel parameters. Thus, the measured ACF of sub-mmWave channel in [29] is chosen. Here, we ignore the effect of beam width and reconfigure some parameters as $f_0$ = 2 GHz, $∥{\mathbf{L}}_{\mathrm{tx}}\left(t_0\right)∥$ = 300 m, K = $7.4$ dB and $∥{\mathbf{v}}_{\mathrm{rx}}\left(t\right)∥$ = $1.2$ m/s according to the measurement condition. Figure 10 gives the comparison between the simulated and measured results and the well agreement reflects the generality and correctness of the proposed model.

5. Conclusions

This paper proposed a 3D non-stationary GBSM for U2V mmWave beam channel by considering the 3D arbitrary trajectory of both terminals, the rotation of 3D antenna array, and 3D beam-forming. To achieve the tradeoff of complexity and accuracy, the model only includes the LoS path and three strongest NLoS paths. Moreover, the hybrid computation method of channel parameters has also been given, which is divided into a geometry-based part and a data-based stochastic part to guarantee both the precision and efficiency. At last, the simulation results of PDPs, ACFs, and DPSDs have been compared with the theoretical and measured ones. In the future, we will perform more channel measurements to verify the channel model as well as optimize the calculation method of channel parameters.