==== Front Sensors (Basel) Sensors (Basel) sensors Sensors (Basel, Switzerland) 1424-8220 MDPI 33291413 10.3390/s20236957 sensors-20-06957 Article A Geometry-Based Beamforming Channel Model for UAV mmWave Communications Mao Kai 1 https://orcid.org/0000-0002-4995-5970Zhu Qiuming 1* https://orcid.org/0000-0001-8183-9139Song Maozhong 1 Hua Boyu 1 Zhong Weizhi 2 Ye Xijuan 1 1 The Key Laboratory of Dynamic Cognitive System of Electromagnetic Spectrum Space, College of Electronic and Information Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, China; maokai@nuaa.edu.cn (K.M.); smz108@nuaa.edu.cn (M.S.); byhua@nuaa.edu.cn (B.H.); yexijuan@nuaa.edu.cn (X.Y.) 2 The Key Laboratory of Dynamic Cognitive System of Electromagnetic Spectrum Space, College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, China; zhongwz@nuaa.edu.cn * Correspondence: zhuqiuming@nuaa.edu.cn 05 12 2020 12 2020 20 23 695727 10 2020 03 12 2020 © 2020 by the authors.2020Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).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. UAVbeamforming channelsGBSMmmWave communicationsray tracing ==== Body 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 (1) vT/R(t)=vT/R(t)cosαT/Rv(t)cosβT/Rv(t)sinαT/Rv(t)cosβT/Rv(t)sinβT/Rv(t) where vT/R(t), αT/Rv(t) and βT/Rv(t) are the magnitude, azimuth angle, and elevation angle of vT/R(t), respectively. It should be noted that ·T/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 hLoS(t) and three strongest NLoS paths, i.e., hGS(τ,t), hSB1(τ,t), hSB2(τ,t) as (2) h(τ,t)=hLoS(t)+hGS(τ,t)︸groundspecular+hSB1(τ,t)+hSB2(τ,t)︸Singlebounce=K(t)K(t)+1A(θLoS,φLoS)h˜LoS(t)+1K(t)+1∑j=13∑m=1MA(θmj,φmj)Pmj(t)h˜mj(t)δ(t−τmj) where M is the ray number of each NLoS path, j∈1,2,3 represents the GS path, SB1 path and SB2 path, respectively. In (2), A(θLoS,φLoS) and A(θmj,φmj) are the gain coefficient of LoS path and NLoS paths within the beam, respectively, Pmj(t) and τmj(t) are the power and delay of the mth ray within the jth NLoS path, respectively. Moreover, h˜mj(t) denotes the channel coefficient of mth ray and can be modeled as (3) h˜mj(t)=exp(jΦI)expj2πrR,mj(t)·RR(t)·dR(t)λ·expj2πrT,mj(t)·RT(t)·dT(t)λ·expj2π∫t0tfmj(t′)dt′ where ΦI represents the random initial phase distributed uniformly over [0,2π), λ is the wavelength, t0 is the initial time instant, and dT/R(t) denotes the location vectors of UAV transmitting antenna (or vehicle receiving antenna) in their own coordinate systems and can be described as (4) dT/R(t)=dT/Rx(t)dT/Ry(t)dT/Rz(t) where dT/Rx(t), dT/Ry(t), and dT/Rz(t) represent the x, y, and z component of dT/R(t). In (3), rT,mj(t) (or rR,mj(t)) and fmj are the spherical unit vectors and Doppler frequency of the mth ray, respectively, and rT,mj(t) (or rR,mj(t)) can be further expressed as (5) rT/R,mj(t)=cosβT/R,mj(t)cosαT/R,mj(t)cosβT/R,mj(t)sinαT/R,mj(t)sinβT/R,mj(t) where αT/R,mj is the azimuth angle of departure (AAoD) or arrival (AAoA), βT/R,mj 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 RT/R(t) is introduced to take this factor into account as (6) RT/R(t)=cosαT/Rv(t)cosβT/Rv(t)−sinαT/Rv(t)−cosαT/Rv(t)sinβT/Rv(t)sinαtx/rxv(t)cosβT/Rv(t)cosαT/Rv(t)−sinαT/Rv(t)sinβT/Rv(t)sinβT/Rv(t)0cosβT/Rv(t) 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 (7) h˜LoS(t)=exp−j2πDLoS(t)λexpj2πrRLoS(t)·RR(t)·dR(t)λ·expj2πrTLoS(t)·RT(t)·dT(t)λexpj2πλ∫0tfLoS(t′))dt′δ(t−τLoS(t)) where DLoS(t) is the distance between the UAV and vehicle, rT/RLoS(t) and fLoS denote the spherical unit vectors and Doppler frequency of LoS path, respectively, and rT/RLoS(t) can be expressed by αT/RLoS and βT/RLoS as (8) rT/RLoS(t)=cosβT/RLoS(t)cosαT/RLoS(t)cosβT/RLoS(t)sinαT/RLoS(t)sinβT/RLoS(t). 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 (9) LT(t)=LT(t0)+∫t0tvT(t)dt (10) LR(t)=LR(t0)+∫t0tvR(t)dt where LT/R(t0) denotes the initial location vector of UAV (or vehicle) at t=t0. The distance vector between the UAV and vehicle DLoS(t) can be expressed as (11) DLoS(t)=LT(t)−LR(t)=DLoS(t0)rTLoS(t)+∫t0tvT,R(t)dt where DLoS(t0) is the initial distance of LoS path and can be expressed as (12) DLoS(t0)=LT(t0)−LR(t0) And vT,R(t) is the relative velocity between the UAV and vehicle which can be expressed as (13) vT,R(t)=vT(t)−vR(t). Similarly, the distance vector between the UAV (or vehicle) and jth scatterer DT/R,j(t) can be expressed as (14) DT/R,j(t)=LT/R(t)−Lj(t)=DT/R,jt0;trT/Rj(t)+∫t0tvT/R(t)dt where DT/R,jt0;t is the initial distance between the UAV (or vehicle) and jth scatterer. In (14), rT/Rj(t) is the spherical unit vectors of each NLoS path, which can be expressed by the mean angles denoted by α¯T/Rj and β¯T/Rj as (15) rT/Rj(t)=cosβ¯T/R,mj(t)cosα¯T/R,mj(t)cosβ¯T/R,mj(t)sinα¯T/R,mj(t)sinβ¯T/R,mj(t) 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 DT/R,jt0;t can be rewritten as (16) DT/R,jt0;t=DT/R,j(t0+kT0)Wt−kT0,k=0,1,2... where the window function W(·) is introduced and can be denoted as (17) Wt=Δ1,0≤t≤T00,otherwise where T0 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 (18) DLoS(t)=DLoS(t0)cos(αT/RLoS(t0))cos(βT/RLoS(t0))+∫t0t(vT,R(t))·cos(αT,Rv(t))·cos(βT,Rv(t))dt2+DLoS(t0)cos(βT/RLoS(t0))sin(αT/RLoS(t0))+∫t0t(vT,R(t))·sin(αT,Rv(t))cos(βT,Rv(t))dt2+DLoS(t0)sin(βT/RLoS(t0))+∫t0t(vT,R(t))·sin(βT,Rv(t))dt2 (19) DT/R,j(t)=DT/R,jt0;tcos(α¯T/Rj(t0))cos(β¯T/Rj(t0))+∫t0tvT/R(t)cos(αT/Rv(t))cos(βT/Rv(t))dt2+DT/R,jt0;tcos(β¯T/Rj(t0))sin(α¯T/Rj(t0))+∫t0tvT/R(t)sin(αT/Rv(t))cos(βT/Rv(t))dt2+DT/R,jt0;tsin(β¯T/Rj(t0))+∫t0tvT/R(t)sin(βT/Rv(t))dt2 where αT/Rv(t) and βT/Rv(t) 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 (20) αT/RLoS(t)=arccos(DT/R,LoSx(t)DT/R,LoSx(t)2+DT/R,LoSy(t)2),DT/R,LoSx(t)≥0π−arccos(DT/R,LoSx(t)DT/R,LoSx(t)2+DT/R,LoSy(t)2),DT/R,LoSx(t)<0 (21) βT/RLoS(t)=arcsin(DT/R,LoSz(t)DLoS(t)) where DT/R,LoSx(t), DT/R,LoSy(t), DT/R,LoSz(t) denote the x, y, z component of DT/R,LoS(t), respectively. Then, the Doppler frequency of LoS path can be rewritten as (22) fLoS(t)=vT(t)·rTLoS(t)λ+vR(t)·rRLoS(t)λ=vT(t)cos(αTLoS(t)−αTv(t))cosβTLoS(t)cosβTv(t)+sinβTLoS(t)sinβTv(t)λ++vR(t)cos(αRLoS(t)−αRv(t))cosβRLoS(t)cosβRv(t)+sinβRLoS(t)sinβRv(t)λ. For the NLoS paths, the mean angles of time-variant EAoD, AAoD or EAoA, AAoA can be calculated respectively by (23) α¯T/Rj(t)=arccos(DT/R,jx(t)DT/R,jx(t)2+DT/R,jy(t)2),DT/R,jx(t)≥0π−arccos(DT/R,jx(t)DT/R,jx(t)2+DT/R,jy(t)2),DT/R,jx(t)<0 (24) β¯T/Rj(t)=arcsin(DT/R,jz(t)DT/R,j(t)). Similarly, the time-variant delays of LoS and NLoS paths are related with the transmission distance and they can be calculated respectively by (25) τLoS(t)=DLoS(t)c (26) τ¯j(t)=DT,j(t)+DR,j(t)c 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 (27) f(Δτmj)=aτ2πστexp(−(Δτmj)22στ2)+bτ where aτ and bτ are the weight factor and correction factor of the magnitude respectively, and στ is the standard deviation of delay offset. In Figure 3, aτ is 9.36, 9.80, 9.56, 10.20, στ is 11.24, 5.10, 8.11, 8.13 and bτ is 0.011, 0.005, 0.015, 0.016 in four different scenarios. Then, the ray delay can be generated by (28) τmj(t)=τ¯j(t)+Δτmj. Furthermore, the power of each ray can be calculated according to the ray delay by the exponential distribution as (29) Pmj(t)=ξexp−1000ξτmj(t) where ξ is the rate parameter of exponential distribution. According to the analytical results of the simulation data, ξ 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 (30) P˜mj(t)=Pmj(t)A(θLoS,φLoS)·K(t)·∑m=1MPmj(t) where the total cluster powers should be 1/(K(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α, σα, and Δαmj as (31) f(Δαmj)=aα2πσαexp(−(Δαmj)22σα2)+bα. In Figure 4a, aα is 0.51, 0.98, 0.68, 0.89, σα is 0.81, 0.51, 0.60, 0.60 and bα is 0.017, 0.006, 0.014, 0.013, respectively. The PDFs of elevation offset EAoA or EAoD can be expressed by Δβmj as (32) f(Δβmj)=aβ2σβexp(−Δβmjσβ). In Figure 4b, aβ is 1.810, 1.008, 1.005, 1.007 and σβ is 3.52, 0.57, 0.70, 0.58 in different scenarios. Then the ray angle can be generated by (33) αT/R,mj(t)=α¯T/Rj(t)+Δαmj (34) βT/R,mj(t)=β¯T/Rj(t)+Δβmj. 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 fmj(t) within jth NLoS path can be described as (35) fmj(t)=fT,mj(t)+fR,mj(t) where fT,mj(t) and fR,mj(t) can be further expressed as (36) fT,mj(t)=vT(t)·rT,mj(t)λ=vT(t)cos(αT,mj(t)−αTv(t))cosβT,mj(t)cosβTv(t)+sinβT,mj(t)sinβTv(t)λ (37) fR,mj(t)=vR(t)·rR,mj(t)λ=vR(t)cos(αR,mj(t)−αRv(t))cosβR,mj(t)cosβRv(t)+sinβR,mj(t)sinβRv(t)λ. 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. 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 t1 = 0 s, t2 = 5 s, and t3 = 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 t1 = 0 s, t2 = 5 s, and t3 = 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 [0∘90∘]. 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 (38) rhhΔt;t=Eh*(t)h(t+Δt)Eh*(t)2Eh(t+Δt)2 where E· is the expectation function and ·* is the complex conjugate. By submitting (2) into (38), 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 t1 = 0 s, t2 = 5 s, and t3 = 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 f0 = 2 GHz, Ltx(t0) = 300 m, K = 7.4 dB and vrx(t) = 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. Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Author Contributions Software and Writing—original draft, K.M.; Conceptualization and funding acquisition, Q.Z.; Methodology, M.S.; Validation, B.H.; Writing—review and editing, W.Z.; Investigation, X.Y. All authors have read and agreed to the published version of the manuscript. Funding This research was funded in part by the National Key Scientific Instrument and Equipment Development Project under Grant No. 61827801, in part by CEMEE State Key Laboratory fund, No. 2020Z0207B, in part by National Defense Science and Technology Key Laboratory fund, No. 6142001190105, in part by Aeronautical Science Foundation of China, No. 201901052001, and in part by the Fundamental Research Funds for the Central Universities, No. NS2020026 and No. NS2020063. Conflicts of Interest The authors declare no conflict of interest. Figure 1 Typical U2V mmWave beam channel. Figure 2 Analyzed results of path power for U2V mmWave channels. Figure 3 The PDFs of ray delay offset. Figure 4 The PDFs of the ray angle offset. Figure 5 The power gain coefficient of 3D beam in different angles. Figure 6 The time-variant normalized PDPs of proposed U2V channel model. Figure 7 The time-variant PAPs of LoS path and NLoS paths. Figure 8 The simulated and theoretical time-variant ACFs at three time instants. Figure 9 The simulated 3D time-variant DPSDs. Figure 10 The simulated and measured results of ACFs. sensors-20-06957-t001_Table 1Table 1 Simulation parameters. Definition Value Definition Value vT(t) 10+0.5t m/s vR(t) 2+t m/s αTv(t) 120−2t∘ βTv(t) 6+5t∘ αRv(t) −120+2t∘ βRv(t) 0∘ λ 3/280 m K 7 dB LT(t0) 400 m LR(t0) 100 m στ 11.24 σα 0.81 σβ 3.52 M 12 ==== Refs References 1. Zeng Y. Wu Q. Zhang R. Accessing from the sky: A tutorial on UAV communications for 5G and beyond Proc. IEEE 2019 107 2327 2375 10.1109/JPROC.2019.2952892 2. Li B. Fei Z. Zhang Y. UAV communications for 5G and beyond: Recent advances and future trends IEEE Internet Things J. 2018 6 2241 2263 10.1109/JIOT.2018.2887086 3. Zhang L. Zhao H. Hou S. Zhao Z. Xu H. Wu X. Wu Q. Zhang R. A survey on 5G millimeter wave communications for UAV-assisted wireless networks IEEE Access 2019 7 117460 117504 10.1109/ACCESS.2019.2929241 4. Zhong W. Xu L. Zhu Q. Chen X. Zhou J. MmWave beamforming for UAV communications with unstable beam pointing China Commun. 2019 16 37 46 10.23919/JCC.2019.10.002 5. Zhong W. Xu L. Liu X. Zhu Q. Zhou J. Adaptive beam design for UAV network with uniform plane array Phys. Commun. 2019 34 58 65 10.1016/j.phycom.2019.02.007 6. Khawaja W. Guvenc I. Matolak D.W. Fiebig U. Schneckenburger N. A survey of air-to-ground propagation channel modeling for unmanned aerial vehicles IEEE Commun. Surv. Tutor. 2020 34 1 6 10.1109/COMST.2019.2915069 7. Cheng X. Li Y. A 3-D geometry-based stochastic model for UAV-MIMO wideband non-stationary channels IEEE Internet Things J. 2018 6 1654 1662 10.1109/JIOT.2018.2874816 8. Cui Z. Briso-Rodriguez C. Guan K. Calvo-Ramirez C. Ai B. Zhong Z. Measurement-based modeling and analysis of UAV air-ground channels at 1 and 4 GHz IEEE Antennas Wirel. Propag. Lett. 2019 18 1804 1808 10.1109/LAWP.2019.2930547 9. Chu X. Briso C. He D. Yin X. Dou J. Channel modeling for low-altitude UAV in suburban environments based on ray tracer Proceedings of the 12th European Conference on Antennas and Propagation (EuCAP 2018) London, UK 9–13 April 2018 1 4 10. Chen X. Hu X. Zhu Q. Zhong W. Chen B. Channel modeling and performance analysis for UAV relay systems China Commun. 2018 15 89 97 11. Goddemeier N. Wietfeld C. Investigation of air-to-air channel characteristics and a UAV specific extension to the Rice model Proceedings of the IEEE Globecom Workshops (GC Wkshps) San Diego, CA, USA 6–10 December 2015 1 5 12. Zheng J. Zhang J. Chen S. Zhao H. Ai B. Wireless powered UAV relay communications over fluctuating two-ray fading channels Phys. Commun. 2019 35 100724 10.1016/j.phycom.2019.100724 13. Chang H. Bian J. Wang C.-X. Aggoune E.M. A 3D non-stationary wideband GBSM for low-altitude UAV-to-ground V2V MIMO channels IEEE Access 2019 7 70719 70732 10.1109/ACCESS.2019.2919790 14. Zhang X. Cheng X. Three-dimensional non-stationary geometry-based stochastic model for UAV-MIMO Ricean fading channels IET Commun. 2018 13 2617 2627 10.1049/iet-com.2019.0139 15. Zhu Q. Jiang K. Chen X. Zhong W. Yang Y. A novel 3D non-stationary UAV-MIMO channel model and its statistical properties China Commun. 2018 15 147 158 16. Zhu Q. Wang Y. Jiang K. Chen X. Zhong W. Ahmed N. 3D non-stationary geometry-based multi-input multi-output channel model for UAV-ground communication systems IET Microw. Antennas Propag. 2019 13 1104 1112 10.1049/iet-map.2018.6129 17. Huang J. Liu Y. Wang C.-X. Sun J. Xiao H. 5G millimeter wave channel sounders, measurements, and models: Recent developments and future challenges IEEE Commun. Mag. 2019 57 138 145 10.1109/MCOM.2018.1701263 18. Wang C.-X. Bian J. Sun J. Zhang W. Zhang M. A survey of 5G channel measurements and models IEEE Commun. Surv. Tutor. 2018 20 3142 3168 10.1109/COMST.2018.2862141 19. Gonzalez-Plaza A. Calvo-Ramirez C. Briso-Rodriguez C. Garcia-Loygorri J.M. Oliva D. Alonso J.I. Propagation at mmW band in metropolitan railway tunnels Wirel. Commun. Mob. Comput. 2017 2018 1 10 10.1155/2018/7350494 20. Liu X. Yin X. Zheng G. Experimental investigation of millimeter-wave MIMO channel characteristics in tunnel IEEE Access 2019 7 108395 108399 10.1109/ACCESS.2019.2932576 21. Bas C.U. Wang R. Sangodoyin S. Hur S. Whang K. Park J. Zhang J. Molisch A.F. 28 GHz propagation channel measurements for 5G microcellular environments Proceedings of the 2018 International Applied Computational Electromagnetics Society Symposium (ACES) Denver, CO, USA 25–29 March 2018 1 2 22. Fu Z. Cui H. Geng S. Zhao X. 5G Millimeter Wave Channel Modeling and Simulations for a High-Voltage Substation Proceedings of the 2019 IEEE Sustainable Power and Energy Conference (iSPEC) Beijing, China 21–23 November 2019 1822 1826 23. Yang G. Zhang Y. He Z. Wen J. Ji Z. Li Y. Machine-learning-based prediction methods for path loss and delay spread in air-to-ground millimetre-wave channels IET Microw. Antennas Propag. 2019 13 1113 1121 10.1049/iet-map.2018.6187 24. Khawaja W. Ozdemir O. Guvenc I. Temporal and spatial characteristics of mm wave propagation channels for UAVs Proceedings of the 2018 11th Global Symposium on Millimeter Waves (GSMM) Boulder, CO, USA 22–24 May 2018 1 6 25. Michailidis E.T. Nomikos N. Trakadas P. Kanatas A.G. Three-dimensional modeling of mmWave doubly massive MIMO aerial fading channels IEEE Trans. Veh. Technol. 2020 69 1190 1202 10.1109/TVT.2019.2956460 26. Cheng L. Zhu Q. Wang C.-X. Zhong W. Hua B. Jiang S. Modeling and simulation for UAV air-to-ground mmWave channels Proceedings of the 2020 14th European Conference on Antennas and Propagation (EuCAP) Copenhagen, Denmark 15–20 March 2020 1 5 27. 3GPP TR 38.901 Study on Channel Model for Frequencies from 0.5 to 100 GHz Available online: ftp://www.3gpp.org/specs/archive/38_series/38.901 (accessed on 17 May 2019) 28. Zhong W. Gu Y. Zhu Q. Li P. Chen X. A novel spatial beam training strategy for mmWave UAV communications Phys. Commun. 2020 41 101106 10.1016/j.phycom.2020.101106 29. Simunek M. Fontan F.P. Pechac P. The UAV low elevation propagation channel in urban areas: Statistical analysis and time-series generator IEEE Trans. Antennas Propag. 2013 61 3850 3858 10.1109/TAP.2013.2256098 30. Cid E.L. Alejos A.V. Sanchez M.G. Signaling through scattered vegetation: Empirical loss modeling for low elevation angle satellite paths obstructed by isolated thin trees IEEE Veh. Technol. Mag. 2016 11 22 28 31. Khawaja W. Ozdemir O. Guvenc I. UAV air-to-ground channel characterization for mmWave systems Proceedings of the 2017 IEEE 86th Vehicular Technology Conference (VTC-Fall) Toronto, ON, Canada 24–27 September 2017 1 5 32. Geise R. Weiss A. Neubauer B. Modulating features of field measurements with a UAV at millimeter wave frequencies Proceedings of the 2018 IEEE Conference on Antenna Measurements & Applications (CAMA) Vasteras, Sweden 3–6 September 2018 1 4