A Non-Stationary Cluster-Based Channel Model for Low-Altitude Unmanned-Aerial-Vehicle-to-Vehicle Communications

Abstract

Under the framework of sixth-generation (6G) wireless communications, the unmanned aerial vehicle (UAV) plays an irreplaceable role in a number of communication systems. In this paper, a novel cluster-based low-altitude UAV-to-vehicle (U2V) non-stationary channel model with uniform planar antenna arrays (UPAs) is proposed. In order to comprehensively model the scattering environment, both single and twin clusters are taken into account. A novel continuous cluster evolution algorithm that integrates time evolution and array evolution is developed to capture channel non-stationarity. In the proposed algorithm, the link between the time evolution of twin clusters and that of single clusters is established to regulate the temporal evolution trend. Moreover, an improved observable radius method is applied to UPAs for the first time to describe array evolution. Based on the combination of cluster evolution and time-variant channel parameters, some vital statistical properties are derived and analyzed, including space–time correlation function (ST-CF), angular power spectrum density (PSD), Doppler PSD, Doppler spread (DS), frequency correlation function (FCF), and delay spread (RS). The non-stationarity in the time, space, and frequency domain is captured. It demonstrates that the airspeed, density of scatterers within clusters, and carrier frequency have an impact on statistical properties. Furthermore, twin clusters have more flexible spatial characteristics with lower power than single clusters. These conclusions can provide assistance and reference for the design and deployment of 6G UAV communication systems.

1. Introduction

The unmanned aerial vehicle (UAV) has been applied to a wide range of fields due to its unique features of flexibility, portability, and strong adaptability. With the development of the sixth-generation (6G) communication network, the structure of the space–air–ground-integrated network (SAGIN) is considered as a key issue [1,2]. The UAV serving as a mobile base station plays an irreplaceable role in building the SAGIN [1,3]. Due to the improved communication efficiency and enhanced transmission performance of the UAV communication, the combination of UAV communications and other key communication technologies can achieve a wide range of applications, including global coverage through the deployment of the SAGIN and Internet of Things [4]. However, although the UAV brings many advantages to communication systems, the research of UAV characteristics is still in its infancy. Unique features that the UAV contains, including highly dynamic communication channel characteristics and obvious spatial and temporal non-stationarity [5], urgently need a comprehensive and in-depth exploration for the design of an efficient UAV communication system. In order to capture and analyze these features, investigating and studying UAV channels is essential.

Currently, research on UAV channels can be classified into channel measurements and channel modelings. Measurement-based methods can provide empirical channel models and UAV channel properties can be investigated based on the measurement data [5]. Most research shows that factors that influence measurement results of UAV channels include the distribution of building sizes, heights, and densities [6]. According to these influencing factors, usual measurement scenarios can be classified into urban, suburban, and rural areas. A measurement is conducted on the campus, and the UAV ascends from an open field area [7]. Large-scale and small-scale channel parameters are extracted from the measurement at 1 GHz and 4 GHz, finding that the shadowing is related to the channel environment and the fast fading presents strong correlations with the carrier frequency. Wang et al. carry out extensive UAV-to-ground measurement campaigns in typical urban macro-cell and rural macro-cell scenarios [8]. The corresponding path loss model with correlation factors is proposed. A measurement campaign is also conducted in a suburban scenario [9] and some measurement results are presented, including path loss, small-scale fading, and power delay profiles. In addition, Matolak et al. make considerable efforts toward the measurement of UAV-to-ground channels in various scenarios [10]. Meanwhile, as a key technology of 6G communications, multiple-input multiple-output (MIMO) has also become increasingly popular. It aims to achieve efficient communication and a large channel capacity. It is demonstrated that the deployment of the MIMO for UAV channels can yield spatial diversity dependent on the surrounding environment [11,12]. MIMO can improve channel capacity markedly [13]. Although accurate channel measurement data and channel models are available using measurement-based methods, they have significant difficulties in obtaining data, including the demand for high-precision equipment, a high time consumption, and the difficulty in measuring in harsh environments.

Different from channel measurement, channel modeling achieves high efficiency by sacrificing accuracy and can adapt to various scenarios. Channel modeling can be further divided into deterministic channel modeling and stochastic channel modeling. The ray-tracing algorithm is an important channel analysis instrument for deterministic channel modeling. Ray-tracing simulations are performed to estimate the received power, time-variant channel impulse response, and delay spread [14]. However, deterministic channel modeling is dependent on a range of environmental parameters, which will increase computation complexity.

In the meanwhile, as a commonly used stochastic modeling method, geometry-based stochastic channel modelings (GBSMs) can greatly increase the computing efficiency due to concise geometric relationships and the high dependence on environmental parameters being greatly reduced. Therefore, more and more research is focusing on the GBSMs of the UAV channel. Compared with prevailing vehicle-to-vehicle (V2V) channels, the modeling of UAV-to-vehicle (U2V) channels must consider more elements, such as the height of the UAV and the movement mode of the UAV in three-dimensional (3D) space. Undoubtedly, it is logical and universal for UAV channels to be modeled as 3D models.

GBSMs can be further subdivided into regular-shaped (RS) GBSMs and cluster-based GBSMs. Thanks to simple geometric characteristics, RS GBSMs are frequently constructed to model various channel environments [15,16,17,18,19,20,21,22,23,24,25,26] and von-Mises distribution is always used to describe the scattering distribution. The height of the UAV, as a unique channel parameter for UAV channels, needs to be considered to ensure the rationality of the modelings. For high-altitude UAV communications, the surrounding area of the UAV is spacious. Therefore, it is illogical to establish a scattering model around the UAV. Bian et al. simulate a high-altitude U2V communication scenario using a multi-cylinder model that covers scatterers in the distance and scatterers around the vehicle [15]. In addition, an aeronautic random model is constructed to analyze the influence of UAV trajectories on channel properties. The hemisphere model is applied to investigate the impact of UAV rotations [16]. A sphere model [17,18] is proposed to simulate the non-isotropic environment around the ground base station. For low-altitude communication scenarios, scatterers around the UAV should be fully taken into consideration. Cylinder models [19,20] are constructed around the UAV and the vehicle. The motion mode of the UAV is a research priority. Other RS GBSMs used to simulate the distribution of scatterers such as ellipsoid models [21,22] and models of multiple geometric combinations [23,24] are also available.

Although RS GBSMs have simply geometric relationships and a high computational efficiency, the distribution of scatterers is concurrently limited by shapes of regular geometries. For complex scattering scenarios, the simulated effect of RS GBSMs is deservedly weak due to the low flexibility of scattering distribution. To remedy the deficiency, cluster-based GBSMs emerge as the times require. Essentially, cluster-based GBSMs can be seen as combinations of multiple independent and identically distributed RS GBSMs [27,28,29,30] and thus can more accurately simulate a complex scattering distribution and provide diverse velocities of scatterers. Due to the notable advantages, the focus is increasingly shifting to cluster-based GBSMs. A certain number of literature has shown that an evolution algorithm is crucial for capturing the channel non-stationarity, and it is generally developed in the time, frequency, and space domain [31,32,33,34,35,36,37,38,39,40,41,42,43]. The evolution algorithm often depends on the birth-and-death (BD) process [31,32,33,34], which can be explained by the fact that clusters may be observable or unobservable for each antenna. Improved algorithms are also developed in the relative literature. Bai et al. develop a space–time (ST) BD algorithm based on the impact of UAV parameters on channels [39] and an improved space–time–frequency (STF) BD algorithm based on the frequency-dependent path gain [40]. Huang et al. incorporate a novel correlated cluster method and K-means clustering algorithm into the BD algorithm [41] and they also propose a new hybrid method combining the parametric method and a visibility region (VR) approach [42].

Although most literature tends to establish RS GBSM for UAV channels [15,16,17,18,19,20,21,22,23,24,25,26], it is difficult for RS GBSMs to accurately describe complex scattering scenarios using specific scattering distributions such as von Mises distribution and cosine distribution. As important features of UAV communication, highly dynamic properties and channel non-stationarity are necessarily captured and researched. However, the non-stationarity is neglected in [16,17,18,19,24,25]. It is proved that when the scatterer escapes the VR of the transmitter (Tx) or the receiver (Rx), the scattering component disappears [44]. Therefore, it is difficult to fully characterize the channel non-stationarity solely by deriving time-varying parameters, as scatterers are always observed [15,21,22,31,45]. In addition, the VR of each antenna is unique for MIMO, which also contributes to the spatial non-stationarity. The effect is not taken into account in [23,24,27]. Moreover, the type of clusters is assumed as unitary [33,34,35,36,37,38,39,40] for cluster-based GBSMs. Clusters are only composed of single clusters for single-bounced (SB) components [31,37] or more frequently assumed twin clusters for multi-bounced (MB) components [33,34,35,36]. However, SB components and MB components cannot be ignored necessarily for urban UAV low-altitude channels due to dense and unevenly distributed buildings and it is difficult for a unique type of cluster to fully characterize the channel propagation properties. In addition, due to the scarcity of research on UAV channels containing multiple types of clusters, there is a limited exploration of time evolution algorithms for different types of clusters and the impact of different types of clusters on UAV channels. Although an UAV channel considering SB and MB components is constructed in [41], it does not compare the channel characteristics of single clusters and twin clusters. The comparison of the characteristics of different clusters is vital for analyzing the contribution of different multipath components to the variation trend of statistical properties. In addition, compared with uniform linear arrays (ULAs), uniform planar arrays (UPAs) can be equipped with more antennas in limited space and greatly enhance space utilization; thus, they have profound research value in the MIMO system [46,47]. Analogously, the scarcity of UPA deployments for UAV channels results in little exploration of array evolution for UPAs. A visual comparison of the latest GBSMs is represented in

Table 1. A comparison of the latest GBSMs.

Models	Channels	Methods for Describing Channel Models	Methods for Capturing Non-Stationarity	Antenna Type
Model in [17,18,19,20,25,26]	UAV communications	Regular geometry	/	uniform linear arrays (ULAs)
Model in [15,21,22]	UAV communications	Regular geometry	Time-variant channel parameters	ULAs
Model in [23,24,27]	UAV communications	Regular geometry	Time-variant channel parameters and time evolution	ULAs
Model in [48]	Vehicle communications	Regular geometry	Time-variant channel parameters and array evolution	ULAs
Model in [32]	Vehicle communications	Single clusters and twin clusters	Time-variant channel parameters and time evolution	ULAs
Model in [31,45]	Vehicle communications	Single clusters	Time-variant channel parameters	ULAs
Model in [35,36]	UAV communications	Twin clusters	Time-variant channel parameters, time and array evolution	ULAs
Model in [46]	Vehicle communications	Regular geometry	Time-variant channel parameters, time and array evolution	uniform planar arrays (UPAs)
Model in [41]	UAV communications	Single clusters and twin clusters	Time-variant channel parameters, time and array evolution	ULAs
The proposed model	UAV communications	Single clusters and twin clusters	Time-variant channel parameters, time and array evolution	UPAs

.

To fill the above gaps, the main contributions are summarized as follows.

(1) A novel cluster-based low-altitude UAV-to-vehicle non-stationary channel model with UPAs is proposed. Single and twin clusters are simultaneously modeled to comprehensively characterize the channel propagation characteristics.

(2) In order to capture the channel non-stationarity, a developed continuous time-array algorithm is proposed. It is the first attempt to establish an evolutionary relationship between single and twin clusters along the time axis. Additionally, we propose an array evolution algorithm for UPAs to satisfy the concentration of antenna sets.

(3) The motion of the UAV, the vehicle, and clusters is considered. Time-variant channel parameters and channel properties are derived. Statistical properties including time–space correlation function (ST-CF), angular power spectrum density (PSD), Doppler PSD, frequency correlation function (FCF), Doppler spread, and delay spread (RS) are simulated and numerically analyzed. The simulated result investigates the impact of the airspeed, concentration parameter, and carrier frequency on channel properties. Additionally, channel characteristics of single clusters and twin clusters are also compared and analyzed.

The rest of the article is structured as follows. Section 2 describes the cluster-based UAV-to-vehicle channel model and time-variant parameters are derived. Section 3 depicts an improved evolution algorithm. Section 4 derives the corresponding statistical properties, simulated results are analyzed in Section 5, and conclusions are drawn in Section 6.

2. A Novel MIMO UAV-to-Vehicle Channel Model

Figure 1 describes the low-altitude U2V channel model framework. The UAV communication scenario in low-altitude urban areas is assumed. In the model, the height of the transmitter (Tx) is $H_{\mathrm{T}}$. The initial vertical distance between the UAV and the receiver (Rx) is D. To accurately represent the scattering components, both SB components and MB components are considered. For single clusters, only a single bounce occurs within the single cluster itself. Additionally, the MB transmission is simplified by a twin cluster, which represents the first bounced cluster and the last bounced cluster. The m-th ray of the $n_{\mathrm{T}}$-th single cluster around the Tx or Rx is defined as $C_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}$ or $C_{{\mathrm{n}}_{\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}$, respectively, which may be scattered by dense buildings or moving vehicles. The n-th single cluster around the Tx/Rx moves with velocity $v_{\mathrm{n}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}}$ and azimuth angle ${\theta }_{\mathrm{n}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}}$. Similarly, the n-th twin cluster moves with velocity $v_{\mathrm{n}}^{\mathrm{MB}}$ and azimuth angle ${\theta }_{\mathrm{n}}^{\mathrm{MB}}$. The velocity of the UAV and the vehicle are defined by $v_{\mathrm{T}}$ and $v_{\mathrm{R}}$, respectively. The azimuth angle and elevation angle of the UAV flight are ${\varphi }_{\mathrm{T}}$ and ${\theta }_{\mathrm{T}}$. Similarly, the azimuth angle of the vehicle is ${\theta }_{\mathrm{R}}$. UPAs with $M_{\mathrm{T}/\mathrm{R}}$ rows and $N_{\mathrm{T}/\mathrm{R}}$ columns are deployed on the Tx/Rx. The transmitting antenna and receiving antenna of the m-th row and n-th column are expressed by $Arra{y}_{(\mathrm{m},\mathrm{n})}^{\mathrm{T}}$ and $Arra{y}_{(\mathrm{m},\mathrm{n})}^{\mathrm{R}}$, respectively. The antenna spacing of the Tx and Rx are ${\delta }_{\mathrm{T}}$ and ${\delta }_{\mathrm{R}}$, respectively.

2.1. UPA Coordinate

The time-variant U2V channel model is defined in a global coordinate system (GCS), which is shown in Figure 1. Considering that it is extremely complex for the UPA to calculate coordinates of antennas in GCS due to its space characteristics, a suitable local coordinate system (LCS) needs to be constructed and the local coordinate of $Arra{y}_{(\mathrm{m},\mathrm{n})}^{\mathrm{T}/\mathrm{R}}$ is expressed as $A_{(\mathrm{m},\mathrm{n})}^{\mathrm{T}/\mathrm{R},\mathrm{LCS}}={(\frac{2m-1}{2}{\delta }_{\mathrm{T}/\mathrm{R}},\frac{2n-1}{2}{\delta }_{\mathrm{T}/\mathrm{R}},0)}^{\mathrm{T}}$. To convert the LCS coordinate to the GCS coordinate, R as the transformation matrix is calculated based on rotation matrices of the x-axis, y-axis, and z-axis, which is expressed as

$$\begin{array}{c}\hfill R=\left(\begin{array}{ccc}cos{\alpha }_{\mathrm{T}/\mathrm{R}}& -sin{\alpha }_{\mathrm{T}/\mathrm{R}}& 0\\ sin{\alpha }_{\mathrm{T}/\mathrm{R}}& cos{\alpha }_{\mathrm{T}/\mathrm{R}}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}cos{\beta }_{\mathrm{T}/\mathrm{R}}& 0& sin{\beta }_{\mathrm{T}/\mathrm{R}}\\ 0& 1& 0\\ -sin{\beta }_{\mathrm{T}/\mathrm{R}}& 0& cos{\beta }_{\mathrm{T}/\mathrm{R}}\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& cos{\gamma }_{\mathrm{T}/\mathrm{R}}& -sin{\gamma }_{\mathrm{T}/\mathrm{R}}\\ 0& sin{\gamma }_{\mathrm{T}/\mathrm{R}}& cos{\gamma }_{\mathrm{T}/\mathrm{R}}\end{array}\right)\end{array}$$

(1)

where ${\alpha }_{\mathrm{T}/\mathrm{R}}$, ${\beta }_{\mathrm{T}/\mathrm{R}}$, and ${\gamma }_{\mathrm{T}/\mathrm{R}}$ denote the rotation angles in the LCS of the Tx/Rx along the x-axis, y-axis, and z-axis. Thus, the converted coordinate vectors of the transmitting UPA and the receiving UPA relative to the origin of the LCS can be derived by $A_{(\mathrm{m},\mathrm{n})}^{\mathrm{T}/\mathrm{R},\mathrm{GCS}}=R{A}_{(\mathrm{m},\mathrm{n})}^{\mathrm{T}/\mathrm{R},\mathrm{LCS}}$.

2.2. Channel Impulse Response

In the proposed model, four components exist, including the line of sight (LoS) component, SB components via clusters generated around the Tx or Rx, and MB components. The channel impulse response (CIR) can be described as

$$\begin{array}{cc}\hfill & h_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}(t,\tau )=h_{\mathrm{pq}}^{\mathrm{LoS}}\left(t\right)\delta (\tau -{\tau }_{\mathrm{LoS}}\left(t\right))+{\displaystyle \sum _{n_{\mathrm{T}}=1}^{N_{\mathrm{T}}\left(t\right)}}{h}_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)\delta (\tau -{\tau }_{n_{\mathrm{T}}}\left(t\right))\hfill \\ \hfill & +{\displaystyle \sum _{n_{\mathrm{R}}=1}^{N_{\mathrm{R}}\left(t\right)}}{h}_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),{\mathrm{n}}_{\mathrm{R}}}^{\mathrm{SB},\mathrm{R}}\left(t\right)\delta (\tau -{\tau }_{n_R}\left(t\right))+{\displaystyle \sum _{n_d=1}^{N_d\left(t\right)}}{h}_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),{\mathrm{n}}_{\mathrm{d}}}^{\mathrm{MB}}\left(t\right)\delta (\tau -{\tau }_{n_d}\left(t\right))\hfill \end{array}$$

(2)

where

$$h_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{LoS}}\left(t\right)=\sqrt{\frac{K}{K+1}}{e}^{-j2\pi d_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}\left(t\right)/\lambda }·e^{j2\pi t{f}_{\mathrm{d},\mathrm{los}}\left(t\right)}$$

(3)

$$\begin{array}{cc}\hfill & h_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),{\mathrm{n}}_{\mathrm{T}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)=\hfill \\ \hfill & \sqrt{\frac{{\eta }_1}{K+1}}\sqrt{P_{{\mathrm{n}}_{\mathrm{T}}}}{\displaystyle \underset{M_{n_T}\to \infty }{lim}}\frac{1}{\sqrt{M_{n_T}}}{\displaystyle \sum _{m_{n_T}=1}^{M_{n_T}}}{e}^{-j2\pi d_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)/\lambda }·e^{j2\pi t{f}_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)+j{\varphi }_1}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & h_{(p,q)(p^{\prime },q^{\prime }),n_R}^{\mathrm{SB},\mathrm{R}}\left(t\right)=\hfill \\ \hfill & \sqrt{\frac{{\eta }_2}{K+1}}\sqrt{P_{{\mathrm{n}}_{\mathrm{R}}}}{\displaystyle \underset{M_{n_R}\to \infty }{lim}}\frac{1}{\sqrt{M_{n_R}}}{\displaystyle \sum _{m_{n_R}=1}^{M_{n_R}}}{e}^{-j2\pi d_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}\left(t\right)/\lambda }·e^{j2\pi t{f}_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}\left(t\right)+j{\varphi }_2}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & h_{(p,q)(p^{\prime },q^{\prime }),n_d}^{\mathrm{MB}}\left(t\right)=\hfill \\ \hfill & \sqrt{\frac{{\eta }_3}{K+1}}\sqrt{P_{{\mathrm{n}}_{\mathrm{d}}}}{\displaystyle \underset{M_{n_d}\to \infty }{lim}}\frac{1}{\sqrt{M_{n_d}}}{\displaystyle \sum _{m_{n_d}=1}^{M_{n_d}}}{e}^{-j2\pi d_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}\left(t\right)/\lambda }·e^{j2\pi t{f}_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}\left(t\right)+j{\varphi }_3}\hfill \end{array}$$

(6)

where K and $\lambda$ denote the Rician factor and wavelength, respectively. $P_{\mathrm{pq},{\mathrm{n}}_{\mathrm{i}}}$ denotes the power of the $n_{\mathrm{i}}$-th cluster. ${\eta }_i$ is the power allocated factor, which satisfies ${\eta }_1+{\eta }_2+{\eta }_3=1$. The initial total number of clusters generated around the Tx and Rx is $N_{\mathrm{T}}\left(0\right)$ and $N_{\mathrm{R}}\left(0\right)$, respectively. The initial number of twin clusters is $N_{\mathrm{d}}\left(0\right)$. The $n_{\mathrm{i}}$-th $(i\in \{s,o,d\left\}\right)$ cluster contains $M_{{\mathrm{n}}_{\mathrm{i}}}$ rays, which approximates to infinity. Moreover, rays within each cluster are assumed time-invariant. ${\tau }_{{\mathrm{n}}_{\mathrm{i}}}\left(t\right)$ denotes the delay of the $n_{\mathrm{i}}$-th cluster and the delays of rays within each cluster are indistinguishable. In addition, $f_{n_{\mathrm{i}},m_{n_{\mathrm{i}}}}^{\mathrm{SB},\mathrm{i}}$ and $d_{n_{\mathrm{i}},m_{n_{\mathrm{i}}}}^{\mathrm{SB},\mathrm{i}}$$(i\in \{T,R\left\}\right)$ denote Doppler frequency and the distance from the $(p,q)$-th transmitting antenna to the $(p^{\prime },q^{\prime })$-th receiving antenna via $C_{{\mathrm{n}}_{\mathrm{i}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{i}}}}^{\mathrm{SB},\mathrm{i}}$. Similarly, $f_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}$ and $d_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}$ denote the Doppler frequency of multi-bounced components and distance from the $(p,q)$-th transmitting antenna to the $(p^{\prime },q^{\prime })$-th receiving antenna through $C_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB}}$.

2.3. Channel Parameters

2.3.1. Initialization Parameters

The velocity vector of the Tx is expressed as $\overrightarrow{v_{\mathrm{T}}}=\left|\right|v_{\mathrm{T}}\left|\right|(cos{\varphi }_{\mathrm{T}}cos{\omega }_{\mathrm{T}},cos{\varphi }_{\mathrm{T}}sin{\omega }_{\mathrm{T}}$, $sin{\varphi }_{\mathrm{T}}{)}^{\mathrm{T}}$.

The velocity vector of the Rx can be written as $\overrightarrow{v_{\mathrm{R}}}=\left|\right|v_{\mathrm{R}}\left|\right|{(cos{\omega }_{\mathrm{R}},sin{\omega }_{\mathrm{R}},0)}^{\mathrm{T}}$.

The cluster’s velocity vector is expressed as $\overrightarrow{v_{\mathrm{n}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}}}=\left|\right|v_{\mathrm{n}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}}\left|\right|(cos{\omega }_{\mathrm{n}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}},$

$sin{\omega }_{\mathrm{n}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}}{,0)}^{\mathrm{T}}$.

The initial position coordinates of the UAV and the vehicle are $\overrightarrow{L_{\mathrm{UAV}}}\left(0\right)={(0,0,H_{\mathrm{T}})}^{\mathrm{T}}$ and $\overrightarrow{L_{\mathrm{vehicle}}}\left(0\right)={(D,0,0)}^{\mathrm{T}}$. The coordinates of the $(p,q)$-th transmitting antenna vector and the $(p^{\prime },q^{\prime })$-th receiving antenna vector can be calculated as

$${\overrightarrow{A}}_{(\mathrm{p},\mathrm{q})}^{\mathrm{T}}\left(0\right)=\left(\begin{array}{c}{x}_{A_{(\mathrm{p},\mathrm{q})}^{\mathrm{T},\mathrm{GCS}}}\\ y_{A_{(\mathrm{p},\mathrm{q})}^{\mathrm{T},\mathrm{GCS}}}\\ z_{A_{(\mathrm{p},\mathrm{q})}^{\mathrm{T},\mathrm{GCS}}}+H_{\mathrm{T}}\end{array}\right){\overrightarrow{A}}_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R}}\left(0\right)=\left(\begin{array}{c}D+x_{A_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R},\mathrm{GCS}}}\\ y_{A_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R},\mathrm{GCS}}}\\ z_{A_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R},\mathrm{GCS}}}\end{array}\right)$$

(7)

where $A_{(\mathrm{p},\mathrm{q})}^{\mathrm{T},\mathrm{GCS}}={(x_{A_{(\mathrm{p},\mathrm{q})}^{\mathrm{T},\mathrm{GCS}}},y_{A_{(\mathrm{p},\mathrm{q})}^{\mathrm{T},\mathrm{GCS}}},z_{A_{(\mathrm{p},\mathrm{q})}^{\mathrm{T},\mathrm{GCS}}})}^{\mathrm{T}}$.

The location vector of the m-th scatterer within the $n_{\mathrm{T}/\mathrm{R}}$-th cluster generated around the Tx or Rx is determined by the distance from the cluster to the Tx/Rx and the angle of departure/arrival (AoD/AoA), which can be written as

$$\begin{array}{c}\hfill {\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(0\right)=\left(\begin{array}{c}{D}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}cos{\gamma }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(0\right)cos{\omega }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(0\right)\\ D_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}cos{\gamma }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(0\right)sin{\omega }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(0\right)\\ D_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}sin{\gamma }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(0\right)\end{array}\right)\end{array}$$

(8)

where ${\gamma }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}$ and ${\omega }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}$ are the elevation angle of departure (EAoD)/the elevation angle of arrival (EAoA) and the azimuth angle of departure (AAoD)/the azimuth angle of arrival (AAoA), respectively. $D_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)$ denotes the distance from the cluster generated around Tx/Rx to the Tx/Rx. Analogously, the EAoD/AAoD of the $n_{\mathrm{T}}$-th cluster generated around the Tx and the EAoA/AAoA of the $n_{\mathrm{R}}$-th cluster generated around the Rx are described as ${\gamma }_{{\mathrm{n}}_{\mathrm{T}}}^{\mathrm{T}}$/${\omega }_{{\mathrm{n}}_{\mathrm{T}}}^{\mathrm{T}}$ and ${\gamma }_{{\mathrm{n}}_{\mathrm{R}}}^{\mathrm{R}}$/${\omega }_{{\mathrm{n}}_{\mathrm{R}}}^{\mathrm{R}}$, respectively, which are randomly generated based on the probability distribution shown in

Table 2. Distribution of relevant channel parameters.

Parameter	Distribution	Mean	Derivation
${\omega }_{{\mathrm{n}}_{\mathrm{T}}}^{\mathrm{T}}$,${\omega }_{{\mathrm{n}}_{\mathrm{R}}}^{\mathrm{R}}$	wrapped Gaussian	$60°,180°$	$20°,30°$
${\gamma }_{{\mathrm{n}}_{\mathrm{T}}}^{\mathrm{T}}$, ${\gamma }_{{\mathrm{n}}_{\mathrm{R}}}^{\mathrm{R}}$	wrapped Gaussian	$-45°,45°$	$30°,30°$
$D_n^T,D_n^R$	exponential	$\frac{1}{30}$,$\frac{1}{30}$	$\frac{1}{900}$,$\frac{1}{900}$

.

Rays are usually assumed as randomly generated following Poisson distribution [37,40,41] and thus channel parameters for each ray can be calculated separately. However, the calculation of channel parameters is extremely complex when extensive rays exist within each cluster. An appropriate probability distribution needs to be applied to describe the distribution of scatterers. Von Mises distribution is widely used to describe the distribution of AoDs in RS GBSMs. A two-dimensional von Mises distribution is applied to imitate both the distribution of elevation angles and azimuth angles, in which elevation angles and azimuth angles are assumed to be mutually independent. The probability density function (PDF) for clusters generated around the Tx is written as

$$\begin{array}{c}\hfill p({\omega }_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{T}},{\gamma }_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{T}})=\frac{exp[k_1 cos({\omega }_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{T}} - {\omega }_{{\mathrm{n}}_{\mathrm{T}}}^{\mathrm{T}}) + k_2 cos({\gamma }_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{T}} - {\gamma }_{{\mathrm{n}}_{\mathrm{T}}}^{\mathrm{T}})]}{{\left(2\pi \right)}^2{I}_0\left(k_1\right)I_0\left(k_2\right)}\end{array}$$

(9)

where $k_1$ and $k_2$ denote concentration parameters, which control the width of probability distribution. Scatterers are more concentrated around the mean angles with the increase in concentration parameters. $I_0(·)$ is the zero-th-order Bessel function. Consequently, the density of scatterers within each cluster can be adjusted by $k_1$ and $k_2$. The PDF for clusters around the Rx shares a similar expression.

2.3.2. Time-Variant Channel Parameters

The time-variant position coordinate of the UAV and the vehicle can be written as ${\overrightarrow{L}}_{\mathrm{UAV}}\left(t\right)={\overrightarrow{L}}_{\mathrm{UAV}}\left(0\right)+\overrightarrow{v_{\mathrm{T}}}t$ and ${\overrightarrow{L}}_{\mathrm{vehicle}}\left(t\right)={\overrightarrow{L}}_{\mathrm{vehicle}}\left(0\right)+\overrightarrow{v_{\mathrm{R}}}t$. The position coordinate of the $n_{\mathrm{T}/\mathrm{R}}$-th cluster can be calculated as ${\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)={\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{T}/\mathrm{R}}\left(0\right)+{\overrightarrow{v}}_n^{\mathrm{SB},\mathrm{T}/\mathrm{R}}$. The position coordinate of the $n_{\mathrm{d}}$-th cluster can be calculated as ${\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}/\mathrm{R}}\left(t\right)={\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}/\mathrm{R}}\left(0\right)+{\overrightarrow{v}}_n^{\mathrm{MB},\mathrm{T}/\mathrm{R}}$. The distance from the $(p,q)$-th transmitting antenna to the $(p^{\prime },q^{\prime })$-th receiving antenna directly can be calculated as $d_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}\left(t\right)={\overrightarrow{A}}_{(\mathrm{p},\mathrm{q})}^{\mathrm{T}}\left(t\right)-{\overrightarrow{A}}_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R}}\left(t\right)$, in which ${\overrightarrow{A}}_{(\mathrm{p},\mathrm{q})}^{\mathrm{T}}\left(t\right)={\overrightarrow{A}}_{(\mathrm{p},\mathrm{q})}^{\mathrm{T}}\left(0\right)+\overrightarrow{v_{\mathrm{T}}}t$ and ${\overrightarrow{A}}_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R}}\left(t\right)={\overrightarrow{A}}_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R}}\left(0\right)+\overrightarrow{v_{\mathrm{R}}}t$. Meanwhile, the distance from the $(p,q)$-th transmitting antenna to the $(p^{\prime },q^{\prime })$-th receiving antenna impinging on the $n_{\mathrm{i}}$-th single cluster around the Tx/Rx can be expressed as $d_{{\mathrm{n}}_{\mathrm{i}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{i}}}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}}\left(t\right)=\left|\right|{\overrightarrow{A}}_{(\mathrm{p},\mathrm{q})}^{\mathrm{T}}\left(t\right)-{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{i}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{i}}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)\left|\right|+\left|\right|{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{i}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{i}}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)-{\overrightarrow{A}}_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R}}\left(t\right)\left|\right|$. Analogously, the distance through the $n_{\mathrm{d}}$-th twin cluster is $d_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB}}\left(t\right)=\left|\right|{\overrightarrow{A}}_{(\mathrm{p},\mathrm{q})}^{\mathrm{T}}\left(t\right)-{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}}\left(t\right)\left|\right|+\left|\right|{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}}\left(t\right)-{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}}\left(t\right)\left|\right|+\left|\right|{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}}\left(t\right)-{\overrightarrow{A}}_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}^{\mathrm{R}}\left(t\right)\left|\right|$.

The Doppler frequency of the LoS component can be expressed as

$$f_{\mathrm{d},\mathrm{LoS}}\left(t\right)=\frac{<{\overrightarrow{L}}_{\mathrm{vehicle}}\left(t\right)-{\overrightarrow{L}}_{\mathrm{UAV}}\left(t\right),\overrightarrow{v_{\mathrm{T}}}-\overrightarrow{v_{\mathrm{R}}}>}{\lambda ∥{\overrightarrow{L}}_{\mathrm{vehicle}}\left(t\right)-{\overrightarrow{L}}_{\mathrm{UAV}}\left(t\right)∥}$$

(10)

The Doppler frequency of the single-bounced component impinging on the m-th scatterer of the $n_{\mathrm{T}}$-th cluster generated around the Tx and the m-th scatterer of the $n_{\mathrm{R}}$-th cluster generated around the Rx can be expressed as

$$\begin{array}{c}\hfill f_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)=\frac{1}{\lambda }(\frac{<{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right),\overrightarrow{v_{\mathrm{T}}}-\overrightarrow{v_{\mathrm{n}}^{\mathrm{T}}}>}{∥{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)∥}+\frac{<{\overrightarrow{L}}_{\mathrm{UAV}}\left(0\right)\mathrm{ }-\mathrm{ }{\overrightarrow{L}}_{\mathrm{vehicle}}\left(t\right)\mathrm{ }+\mathrm{ }{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right),\overrightarrow{v_{\mathrm{R}}}\mathrm{ }-\mathrm{ }\overrightarrow{v_{\mathrm{n}}^{\mathrm{T}}}>}{∥{\overrightarrow{L}}_{\mathrm{UAV}}\left(0\right)\mathrm{ }-\mathrm{ }{\overrightarrow{L}}_{\mathrm{vehicle}}\left(t\right)\mathrm{ }+\mathrm{ }{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)∥})\end{array}$$

(11)

$$\begin{array}{c}\hfill f_{{\mathrm{n}}_{\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}\left(t\right)=\frac{1}{\lambda }(\frac{<{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}\left(t\right),\overrightarrow{v_{\mathrm{R}}}\mathrm{ }-\mathrm{ }\overrightarrow{v_{\mathrm{n}}^{\mathrm{R}}}>}{∥{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}\left(t\right)∥}+\frac{<{\overrightarrow{L}}_{\mathrm{UAV}}\left(t\right)\mathrm{ }-\mathrm{ }{\overrightarrow{L}}_{\mathrm{vehicle}}\left(0\right)\mathrm{ }+\mathrm{ }{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}\left(t\right),\overrightarrow{v_{\mathrm{R}}}\mathrm{ }-\mathrm{ }\overrightarrow{v_{\mathrm{n}}^{\mathrm{R}}}>}{∥{\overrightarrow{L}}_{\mathrm{UAV}}\left(t\right)\mathrm{ }-\mathrm{ }{\overrightarrow{L}}_{\mathrm{vehicle}}\left(0\right)\mathrm{ }+\mathrm{ }{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}\left(t\right)∥})\end{array}$$

(12)

The Doppler frequency of the $m_{{\mathrm{n}}_{\mathrm{d}}}$-th double-bounced component can be expressed as

$$\begin{array}{c}\hfill f_{\mathrm{pq},{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB}}\left(t\right)=\frac{1}{\lambda }(\frac{<{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}}\left(t\right),\overrightarrow{v_{\mathrm{T}}}\mathrm{ }-\mathrm{ }\overrightarrow{v_{\mathrm{n}}^{\mathrm{T}}}>}{∥{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}}\left(t\right)∥}+\frac{<{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}}\left(t\right),\overrightarrow{v_{\mathrm{R}}}\mathrm{ }-\mathrm{ }\overrightarrow{v_{\mathrm{n}}^{\mathrm{R}}}>}{∥{\overrightarrow{C}}_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}}\left(t\right)∥})\end{array}$$

(13)

where $<·,·>$ denotes the inner product operator. Moreover, different from single-bounced components, a virtual link [49] exists between the first cluster and the last cluster in MB components. The delay of the $n_{\mathrm{d}}$-th twin cluster can be computed as ${\tau }_{\mathrm{pq},{\mathrm{n}}_{\mathrm{d}}}^{\mathrm{db}}\left(t\right)=d_{{\mathrm{n}}_{\mathrm{d}}}^{\mathrm{MB}}\left(t\right)/c+\widehat{{\tau }_{\mathrm{link}}}$, in which $\widehat{{\tau }_{\mathrm{link}}}$ denotes the virtual delay and is randomly drawn following $\widehat{{\tau }_{\mathrm{link}}}=-r_{\tau }{\delta }_{\tau }ln{u}_n$. $r_{\tau }$ is the delay scalar, ${\delta }_{\tau }$ is a stochastic delay spread, and $u_{\mathrm{n}}$ follows uniform distribution, i.e., $u_{\mathrm{n}}\sim \mathrm{U}(0,1)$. $d_{{\mathrm{n}}_{\mathrm{d}}}^{\mathrm{MB}}\left(t\right)$ denotes the distance from the Tx to the Rx via the $n_{\mathrm{d}}$-th twin cluster. Although clusters only have two states, observable and unobservable, the change in power is a continuous process. Therefore, traditional methods used to simulate the variation in power are obviously unreasonable [50]. In order to describe the smooth transition process when clusters are born or die, the transition region is adopted as follows [51]:

$$\begin{array}{c}\hfill \xi \left(t\right)=\frac{1}{2}-\frac{1}{\pi }arctan(\frac{2\left[L_c+(\left|2t-T_n\right|-T_n)\left|\right|\overrightarrow{v_T}-\overrightarrow{v_R}\left|\right|\right]}{\sqrt{\lambda L_c}})\end{array}$$

(14)

where $T_n$ is the lifetime of the n-th cluster. $L_c$ is the duration of the transition process.

The power of each cluster is jointly affected by the channel environment and the position of the cluster. The power of the n-th cluster can be calculated as

$$\begin{array}{c}\hfill {\widehat{P}}_n\left(t\right)=\xi \left(t\right)e^{-{\tau }_n^i\left(t\right)\frac{r_{\tau }-1}{r_{\tau }{\sigma }_{\tau }}}{10}^{-\frac{X}{10}}\end{array}$$

(15)

where X is the shadowing term following Gaussian distribution, i.e., $X\sim N(0,1)$. i denotes the type of the cluster. The normalized power of the n-th cluster is obtained by

$$\begin{array}{c}\hfill P_n\left(t\right)=\frac{{\widehat{P}}_n\left(t\right)}{{\sum }_{n=1}^N{\widehat{P}}_n\left(t\right)}\end{array}$$

(16)

where N denotes the sum of single clusters and twin clusters.

3. Cluster Evolution

The space–time-frequency (STF) non-stationarity is regarded as a crucial channel characteristic. Currently, the method of capturing non-stationary channels remains an open topic. As shown in Table 1, RS GBSMs always capture the non-stationarity by the motion of the Tx, the Rx, and the regular-shaped geometries [15,21,22]. In these models, it is obviously unreasonable that scatterers are always observable, no matter how long it takes. To overcome this deficiency, a two-state Markov process [24] and a BD method [23] are introduced to capture the dynamic characteristic. However, RS GBSMs undoubtedly have rigorous limitations on the distribution and velocities of scatterers. To address the above issues, a cluster-based channel model is often accompanied by a cluster evolution algorithm to capture the non-stationarity. The combination of time-variant channel parameters and cluster evolution algorithm can largely exhibit channel non-stationarity. To ensure a coherent evolutionary process for clusters of different types and states, time evolution and array evolution are presumed as two consecutive and related processes in the proposed algorithm. Furthermore, due to the motion of the UAV, the vehicle, and clusters, as well as the deployment of UPAs, the developed evolution algorithm is carried out in the time and space domain. Considering that new clusters scattered by SB components or MB components are generated continuously over time, the type of new cluster needs to be distinguished. Therefore, the evolution process of cluster types, which is essentially the trend of single- and twin-cluster generation, is also included in the time evolution of the developed algorithm. In addition, array evolution can be simply described as different cluster sets for different antennas. Similar to the array evolution applied to ULAs, due to the VR of each antenna [42], the array evolution in UPAs also manifests as continuity. The flowchart of time-array evolution is shown in Figure 2.

3.1. Time Evolution and Type Evolution

In the algorithm, a general BD process is introduced to describe the evolution on the time axis for surviving single clusters and twin clusters. The initial number of single clusters and twin clusters is $N_{\mathrm{SB}}\left(0\right)$ and $N_{\mathrm{MB}}\left(0\right)$, respectively. ${\lambda }_{\mathrm{B}}$ and ${\lambda }_{\mathrm{D}}$, representing the average birth rate and death rate of clusters, respectively, always associate with the velocity of the Tx and Rx [40] and the channel environment. The survival probabilities of single clusters and twin clusters obey exponential distribution, which can be written as $P_{\mathrm{sur}}^{\mathrm{T}}(\Delta t)=e^{-\frac{{\lambda }_{\mathrm{D}}(v_T+v_R+P_{\mathrm{T}}{v}_{\mathrm{n}}^{\mathrm{T}})\Delta t}{D_{\mathrm{c}1}}}$, $P_{\mathrm{sur}}^{\mathrm{R}}(\Delta t)=e^{-\frac{{\lambda }_{\mathrm{D}}(v_{\mathrm{T}}+v_{\mathrm{R}}+P_{\mathrm{R}}{v}_{\mathrm{n}}^{\mathrm{R}})\Delta t}{D_{\mathrm{c}2}}}$ and $P_{\mathrm{sur}}^{\mathrm{MB}}(\Delta t)=e^{-\frac{{\lambda }_{\mathrm{D}}(v_{\mathrm{T}}+v_{\mathrm{R}}+P_{\mathrm{MB}}(v_{\mathrm{n}}^{\mathrm{MB},\mathrm{T}}+v_{\mathrm{n}}^{\mathrm{MB},\mathrm{R}}))\Delta t}{D_{\mathrm{c}3}}}$, in which $P_{\mathrm{i}}(i\in \{T,R,MB\})$ represents the proportion of moving clusters and $D_i$ is a scenario-dependent factor. New single clusters are generated obeying Poisson distribution, which is written as $E\left(N_{\mathrm{i}}(t+\Delta t)\right)=\frac{{\lambda }_B}{{\lambda }_D}(1-P_{\mathrm{sur}}^{\mathrm{i}})(i\in \{T,R\})$ [39]. The initial number of single clusters around the Tx or Rx is calculated as $N\left(t_0\right)=\frac{{\lambda }_{\mathrm{B}}}{{\lambda }_{\mathrm{D}}}$. In urban communication scenarios with dense buildings, a large number of multipath components will be concretized into single clusters and twin clusters. As clusters evolve on the time axis, newly generated clusters may contribute to MB components or SB components. Since MB components can be regarded as transmissions through twin clusters, the combination of two newly generated single clusters is feasible to describe a new MB component. Therefore, when new clusters are generated both around the Tx and Rx simultaneously, they have a probability of being related and contributing to the MB component. Additionally, the first cluster and the last cluster in a new MB component correspond to a new cluster generated around the Tx and a new cluster generated around the Rx, respectively. As a consequence, if both clusters around the Tx and Rx produce simultaneously in an updated interval $\Delta t$, they are randomly paired and emerge into twin clusters based on the correlation probability $P_{twin}$.

According to the time evolution process, new clusters are constantly emerging with time. For a newly generated single cluster $C_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}}^{\mathrm{SB},\mathrm{T}/\mathrm{R}}$, the location vector is determined based on the distribution listed in Table 2. And, it can be expressed as

$$\begin{array}{cc}\hfill & {\overrightarrow{C}}_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)=\left(\begin{array}{c}{D}_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)cos{\gamma }_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)cos{\omega }_{\mathrm{n}+1}^{\mathrm{T}/\mathrm{R}}\left(t\right)\\ D_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)cos{\gamma }_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)sin{\omega }_{{\mathrm{n}+1}_{\mathrm{T}}/\mathrm{R}}^{\mathrm{T}/\mathrm{R}}\left(t\right)\\ D_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)sin{\gamma }_{{\mathrm{n}+1}_{\mathrm{T}/\mathrm{R}}}^{\mathrm{T}/\mathrm{R}}\left(t\right)\end{array}\right)\hfill \end{array}$$

(17)

Meanwhile, the location coordinates of survival clusters update uninterruptedly based on geometric relationships in each update time interval $\Delta t$, which are derived in Section 2.3.2.

The birth and death of clusters as well as time-variant channel parameters can cause dynamic power variation. An example of the power variation of clusters is shown in Figure 3. It can be observed that with the introduction of the transition region, the BD property is always accompanied by a smooth increase and decrease in the power.

3.2. Array Evolution

In principle, each antenna has a unique VR of similar size and it has been proven that as the Tx and Rx move into the VR of a cluster concurrently, the cluster can be observed. However, each antenna of the MIMO should be judged due to spatial differences. There is no doubt that antennas farther apart share a smaller VR [43]. As a result, a number of clusters are not always visible to the whole MIMO, which leads to a unique cluster set for each antenna. Correspondingly, closer antennas share a larger VR, which results in clusters observed by limited space continuous antennas. Homoplastically, UPAs that can observe a specific cluster are located within a certain area to satisfy the continuity. To capture the channel non-stationarity in the space domain, the BD property can also be used to describe array evolution for UPAs. Nevertheless, observable antennas with continuity and regionality should also be considered. To simulate the behavior of UPA evolution, the observable radius is applied to the UPAs and the developed algorithm fully considers the time-varying process of antenna sets for survival clusters and new clusters.

For new clusters, the algorithm is performed as follows.

Step 1: Select initial transmitting antenna index $(p,q)$ and receiving antenna index $(p^{\prime },q^{\prime })$ obeying uniform distribution, i.e., $\left[p\right]\sim \mathrm{U}(1,M_T),\left[q\right]\sim \mathrm{U}(1,N_T),\left[p^{\prime }\right]\sim \mathrm{U}(1,M_R),\left[q^{\prime }\right]\sim \mathrm{U}(1,N_R)$, in which $[·]$ denotes the integer function.

Step 2: Calculate the observable radius for each cluster based on exponential distribution, i.e., $R_{\mathrm{T}/\mathrm{R}}\sim \mathrm{E}(e^{\frac{{\lambda }_d}{D_{a1/a2}}})$, in which $D_{a1/a2}$ denotes an array-dependent factor.

Step 3: Calculate the distance from each transmitting antenna index $(x_T,y_T)$ to the initial selected transmitting antenna $(p,q)$ and from each receiving antenna index $(x_R,y_R)$ to the initial selected receiving antenna $(p,q)$. The distance can be derived as follows:

$$\begin{array}{cc}\hfill & d_{(\mathrm{p},\mathrm{q})({\mathrm{x}}_{\mathrm{T}},{\mathrm{y}}_{\mathrm{T}})}^{\mathrm{T}}=\sqrt{{(x_T-p)}^2+{(y_T-q)}^2}{\delta }_{\mathrm{T}}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & d_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })({\mathrm{x}}_{\mathrm{R}},{\mathrm{y}}_{\mathrm{R}})}^{\mathrm{T}}=\sqrt{{(x_R-p^{\prime })}^2+{(y_R-q^{\prime })}^2}{\delta }_{\mathrm{R}}\hfill \end{array}$$

(19)

Step 4: Compare values. If $d_{(\mathrm{p},\mathrm{q})({\mathrm{x}}_{\mathrm{T}},{\mathrm{y}}_{\mathrm{T}})}^{\mathrm{T}}\le R_{\mathrm{T}}$ and $d_{({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })({\mathrm{x}}_{\mathrm{R}},{\mathrm{y}}_{\mathrm{R}})}^{\mathrm{R}}\le R_{\mathrm{R}}$, the cluster can be observed by the $(x_T,y_T)$-th transmitting antenna and the $(x_R,y_R)$-th receiving antenna simultaneously.

Due to the interconnected nature of time evolution and array evolution, the continuity of time evolution can also lead to the continuity of array evolution. The gradient array evolution process can effectively avoid the occurrence of a completely different antenna set for survival clusters within an updated time $\Delta t$. However, the time-varying process of antenna sets for survival clusters is neglected in [39]. In order to simulate the gradient process of antenna sets, the initial antenna index at time point $t+\Delta t$ is randomly selected from the antenna set observing the cluster at time point t based on the uniform distribution. Subsequently, Step 2 and Step 3 are carried out to hunt for the antenna sets at time point $t+\Delta t$.

According to the array algorithm, Figure 4 shows that the array evolution can effectively simulate the distribution of antenna sets observing different clusters. It can be seen that each cluster has a unique antenna set and these antennas exhibit regional distribution. Through these means, the space non-stationarity for UPAs is captured sufficiently.

4. Channel Statistical Properties

In this section, the corresponding statistical properties are deduced.

4.1. Space–Time Correlation Function

The ST-CF of the link between the $(p,q)$-th antenna and the $(p^{\prime },q^{\prime })$-th antenna as well as the link between the $(m,n)$-th antenna and the $(m^{\prime },n^{\prime })$-th antenna is computed by $R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}({\delta }_{\mathrm{T}},{\delta }_{\mathrm{R}},t,\Delta t)=\frac{E\left[h_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}\left(t\right)h_{(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}^{\ast }(t+\Delta t)\right]}{\sqrt{E\left[{\left|h_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}\left(t\right)\right|}^2\right]E\left[{\left|h_{(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t+\Delta t)\right|}^2\right]}}$, in which ${(·)}^{\ast }$ denotes the complex conjugate operation. Due to the introduction of the BD process, the survival probability of clusters becomes smaller as the update time progresses. And, the disappearance of clusters naturally leads to the disappearance of the corresponding channel characteristics. Thus, $P_{\mathrm{sur}}^{\mathrm{i}}$ is a crucial part that contributes to the time correlation to describe the probability of clusters surviving in a time interval. Moreover, the ST-CF can be decomposed into the ST-CFs for different subcomponents, which is due to the mutual independence. The ST-CFs can be deduced as (20)–(23).

$$\begin{array}{cc}\hfill & R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}^{\mathrm{LoS}}({\delta }_{\mathrm{T}},{\delta }_{\mathrm{R}},t,\Delta t)\hfill \\ \hfill & =\frac{K}{K+1}{e}^{-j\frac{2\pi }{\lambda }(d_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}\left(t\right)-d_{(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t+\Delta t))}{e}^{j2\pi t(f_{\mathrm{d},\mathrm{LoS}}\left(t\right)-f_{\mathrm{d},\mathrm{LoS}}(t+\Delta t))-j2\pi f_{\mathrm{d},\mathrm{LoS}}(t+\Delta t)\Delta t}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}^{\mathrm{SB},\mathrm{T}}({\delta }_{\mathrm{T}},{\delta }_{\mathrm{R}},t,\Delta t)\hfill \\ \hfill & =\frac{P_{\mathrm{sur}}^{\mathrm{T}}(\Delta t){\eta }_1}{K+1}{\displaystyle \sum _{n_T=1}^{N_T(t)\cap N_T(t+\Delta t)}}{\displaystyle \underset{M_{n_T}\to \infty }{lim}}\frac{P_{{\mathrm{n}}_{\mathrm{T}}}}{M_{n_T}}E({\displaystyle \sum _{m_{n_T}=1}^{M_{n_T}}}{e}^{-j2\pi (d_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}t-d_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t+\Delta t))/\lambda }\hfill \\ \hfill & e^{j2\pi t(f_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t)-f_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t+\Delta t))}{e}^{j2\pi \Delta t{f}_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t+\Delta t)})\hfill \\ \hfill & =\frac{P_{\mathrm{sur}}^{\mathrm{T}}(\Delta t){\eta }_1{P}_{{\mathrm{n}}_{\mathrm{T}}}}{K+1}{\displaystyle \sum _{n_T=1}^{N_T(t)\cap N_T(t+\Delta t)}}{\displaystyle {\int }_{-\pi }^{\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}{e}^{-j2\pi (d_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}t-d_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t+\Delta t))/\lambda }{e}^{j2\pi t(f_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t)-f_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t+\Delta t))}\hfill \\ \hfill & e^{j2\pi \Delta t{f}_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{SB},\mathrm{T}}(t+\Delta t)}p({\omega }_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{T}},{\gamma }_{{\mathrm{n}}_{\mathrm{T}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{T}}}}^{\mathrm{T}})d({\omega }_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{T}})d({\gamma }_{n_{\mathrm{T}},m_{n_{\mathrm{T}}}}^{\mathrm{T}})\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}^{\mathrm{SB},\mathrm{R}}({\delta }_{\mathrm{T}},{\delta }_{\mathrm{R}},t,\Delta t)\hfill \\ \hfill & =\frac{P_{\mathrm{sur}}^{\mathrm{R}}(\Delta t){\eta }_2{P}_{{\mathrm{n}}_{\mathrm{R}}}}{K+1}{\displaystyle \sum _{n_R=1}^{N_R(t)\cap N_R(t+\Delta t)}}{\displaystyle {\int }_{-\pi }^{\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}{e}^{-j2\pi (d_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}t-d_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}(t+\Delta t))/\lambda }{e}^{j2\pi t(f_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}(t)-f_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}(t+\Delta t))}\hfill \\ \hfill & e^{j2\pi \Delta t{f}_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{SB},\mathrm{R}}(t+\Delta t)}p({\omega }_{{\mathrm{n}}_{\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{R}}}}^{\mathrm{R}},{\gamma }_{{\mathrm{n}}_{\mathrm{R}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{R}}}}^{\mathrm{R}})d({\omega }_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{R}})d({\gamma }_{n_{\mathrm{R}},m_{n_{\mathrm{R}}}}^{\mathrm{R}})\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}^{\mathrm{MB}}({\delta }_{\mathrm{T}},{\delta }_{\mathrm{R}},t,\Delta t)\hfill \\ \hfill & =\frac{P_{\mathrm{sur}}^{\mathrm{MB}}(\Delta t){\eta }_3{P}_{{\mathrm{n}}_{\mathrm{MB}}}}{K+1}{\displaystyle \sum _{n_d=1}^{N_d(t)\cap N_d(t+\Delta t)}}{\displaystyle {\int }_{-\pi }^{\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}{e}^{-j2\pi (d_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}t-d_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}(t+\Delta t))/\lambda }\hfill \\ \hfill & e^{j2\pi t(f_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}(t)-f_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}(t+\Delta t))}{e}^{j2\pi \Delta t{f}_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB}}(t+\Delta t)}p({\omega }_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}},{\gamma }_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}})p({\omega }_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}},{\gamma }_{{\mathrm{n}}_{\mathrm{d}},{\mathrm{m}}_{{\mathrm{n}}_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}})\hfill \\ \hfill & d({\omega }_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}})d({\gamma }_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{T}})d({\omega }_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}})d({\gamma }_{n_{\mathrm{d}},m_{n_{\mathrm{d}}}}^{\mathrm{MB},\mathrm{R}})\hfill \end{array}$$

(23)

4.2. Angular Power Spectrum

The distribution of the power in the angular domain is reflected by the angular power spectrum. In the model, the angular power of the AoDs of clusters is derived, and it can be expressed as

$$\begin{array}{cc}\hfill & \lambda (t,\mu )=\frac{K}{K+1}\delta (\mu -{\mu }_{\mathrm{LoS}}\left(t\right))+\frac{{\eta }_1}{K+1}{\displaystyle \sum _{n_T=1}^{N_T\left(t\right)}}{P}_{{\mathrm{n}}_{\mathrm{T}}}\left(t\right)\delta (\mu -{\mu }_{{\mathrm{n}}_{\mathrm{T}}}\left(t\right))\hfill \\ \hfill & +\frac{{\eta }_2}{K+1}{\displaystyle \sum _{n_R=1}^{N_R\left(t\right)}}{P}_{{\mathrm{n}}_{\mathrm{R}}}\left(t\right)\delta (\mu -{\mu }_{{\mathrm{n}}_{\mathrm{R}}}\left(t\right))+\frac{{\eta }_3}{K+1}{\displaystyle \sum _{n_d=1}^{N_d\left(t\right)}}{P}_{{\mathrm{n}}_{\mathrm{d}}}\left(t\right)\delta (\mu -{\mu }_{{\mathrm{n}}_{\mathrm{d}}}\left(t\right))\hfill \end{array}$$

(24)

where ${\mu }_{{\mathrm{n}}_{\mathrm{T}/\mathrm{R}/\mathrm{d}}}$ denotes the AoD of the corresponding cluster.

4.3. Root Mean Square Delay Spread

RMS-DS is one of the most important evaluation indices used to calculate the delay dispersion, and can be derived as

$$\begin{array}{c}\hfill {\delta }_{\tau }=\sqrt{{\displaystyle \sum _{\mathrm{l}=1}^{\mathrm{L}}}{P}_{\mathrm{l}}\left(t\right){\tau }_{\mathrm{l}}^2-{{\displaystyle \sum _{\mathrm{l}=1}^{\mathrm{L}}}{P}_{\mathrm{l}}\left(t\right){\tau }_{\mathrm{l}}}^2}\end{array}$$

(25)

where $P_{\mathrm{l}}\left(t\right)$ and ${\tau }_{\mathrm{l}}$ are the normalized power and delay of the l-th cluster.

4.4. Doppler Spread

The instantaneous frequency can provide a power distribution in the frequency domain. The instantaneous mean Doppler shift and the instantaneous mean Doppler spread can be calculated by [52]

$$\begin{array}{c}\hfill B_{\mathrm{f}}^{\left(1\right)}=\frac{K}{K+1}{f}_{\mathrm{d},\mathrm{LoS}}+{\displaystyle \sum _{\mathrm{i}=1}^3}{\displaystyle \sum _{\mathrm{n}=1}^{\mathrm{N}}}{\displaystyle \sum _{{\mathrm{m}}_{\mathrm{n}}=1}^{{\mathrm{M}}_{\mathrm{n}}}}\frac{{\eta }_i}{K+1}{P}_{{\mathrm{m}}_{\mathrm{n}}}\left(t\right)f_{\mathrm{n},{\mathrm{m}}_{\mathrm{n}}}\end{array}$$

(26)

$$\begin{array}{c}\hfill B_{\mathrm{f}}^{(2)}=\sqrt{(\frac{K}{K+1}{f}_{\mathrm{d},\mathrm{LoS}}^2+{\displaystyle \sum _{\mathrm{i}=1}^3}{\displaystyle \sum _{\mathrm{n}=1}^{\mathrm{N}}}{\displaystyle \sum _{{\mathrm{m}}_{\mathrm{n}}=1}^{{\mathrm{M}}_{\mathrm{n}}}}\frac{{\eta }_i}{K+1}{P}_{{\mathrm{m}}_{\mathrm{n}}}(t)f_{\mathrm{n},{\mathrm{m}}_{\mathrm{n}}}^2)-{(B_{\mathrm{f}}^{(1)})}^2}\end{array}$$

(27)

4.5. Doppler PSD

Doppler PSD can be calculated by the Fourier transform of $R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t,\Delta t)$ on $\Delta t$, which can be written as

$$\begin{array}{c}\hfill S_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t,f_{\mathrm{D}})={\displaystyle {\int }_{-\infty }^{\infty }}{R}_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t,\Delta t)e^{-j2\pi f_{\mathrm{D}}\Delta t}d\Delta t\end{array}$$

(28)

where $f_{\mathrm{D}}$ denotes the Doppler frequency.

4.6. Frequency Correlation Function

The time-variant channel transfer function $H_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t,f)$ can be derived by the Fourier transform of the CIR. The FCF can be acquired by

$$\begin{array}{c}\hfill R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t,\Delta t,f,\Delta f)=\frac{E[H_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}(t,f)H_{(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}^{\ast }(t+\Delta t,f+\Delta f)]}{\sqrt{E[{\left|H_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime })}(t,f)\right|}^2]E[{\left|H_{(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t+\Delta t,f+\Delta f)\right|}^2]}}\end{array}$$

(29)

The FCF can be derived as

$$\begin{array}{cc}\hfill & R_{(\mathrm{p},\mathrm{q})({\mathrm{p}}^{\prime },{\mathrm{q}}^{\prime }),(\mathrm{m},\mathrm{n})({\mathrm{m}}^{\prime },{\mathrm{n}}^{\prime })}(t,\Delta t,f,\Delta f)=\frac{K}{K+1}{e}^{j2\pi f({\tau }_{\mathrm{LoS}}(\Delta t+t)-{\tau }_{\mathrm{LoS}}\left(t\right))+\Delta f{\tau }_{\mathrm{LoS}}(\Delta t+t)}\hfill \\ \hfill & +\frac{{\eta }_1{P}_{{\mathrm{n}}_{\mathrm{T}}}}{K+1}{\displaystyle \sum _{n_T=1}^{N_T}}{e}^{j2\pi f({\tau }_{{\mathrm{n}}_{\mathrm{T}}}(\Delta t+t)-{\tau }_{{\mathrm{n}}_{\mathrm{T}}}\left(t\right))+\Delta f{\tau }_{{\mathrm{n}}_{\mathrm{T}}}(\Delta t+t)}+\frac{{\eta }_2{P}_{{\mathrm{n}}_{\mathrm{R}}}}{K+1}{\displaystyle \sum _{n_R=1}^{N_R}}{e}^{j2\pi f({\tau }_{{\mathrm{n}}_{\mathrm{R}}}(\Delta t+t)-{\tau }_{{\mathrm{n}}_{\mathrm{R}}}\left(t\right))+\Delta f{\tau }_{{\mathrm{n}}_{\mathrm{R}}}(\Delta t+t)}\hfill \\ \hfill & +\frac{{\eta }_3{P}_{{\mathrm{n}}_{\mathrm{d}}}}{K+1}{\displaystyle \sum _{n_d=1}^{N_d}}{e}^{j2\pi f({\tau }_{{\mathrm{n}}_{\mathrm{d}}}(\Delta t+t)-{\tau }_{{\mathrm{n}}_{\mathrm{d}}}\left(t\right))+\Delta f{\tau }_{{\mathrm{n}}_{\mathrm{d}}}(\Delta t+t)}\hfill \end{array}$$

(30)

5. Result and Analysis

In this section, the statistical properties are simulated and numerically analyzed. The channel parameters are set based on a reasonable assumption. Some of the parameters are given in Table 2, and other channel parameters are set as follows: $f=3.5$ GHz, $H_{\mathrm{T}}=30$ m, $D=300$ m, $v_{\mathrm{n}}^{\mathrm{T}/\mathrm{R}/\mathrm{MB}}\sim \mathrm{U}(0,10)$ m/s, ${\theta }_{\mathrm{n}}^{\mathrm{T}/\mathrm{R}/\mathrm{MB}}\sim \mathrm{U}(0,\pi )$, $v_{\mathrm{T}}=20$ m/s, $v_{\mathrm{R}}=20$ m/s, ${\varphi }_{\mathrm{T}}=\frac{\pi }{4}$, ${\theta }_{\mathrm{T}}=0$, ${\theta }_{\mathrm{R}}=\pi$, ${\delta }_{\mathrm{T}}=\frac{\lambda }{2}$, ${\delta }_{\mathrm{R}}=\frac{\lambda }{2}$, ${\alpha }_{T/R}=\frac{\pi }{4}$, ${\beta }_{\mathrm{T}/\mathrm{R}}=\frac{\pi }{4}$, ${\gamma }_{\mathrm{T}/\mathrm{R}}=\frac{\pi }{3}$, $r_{\tau }=2.3$, ${\delta }_{\tau }={10}^{-6.63}$ [39], $M_{\mathrm{T}/\mathrm{R}}=16$, $N_{\mathrm{T}/\mathrm{R}}=16$. For low-altitude UAV communications, the environments around the UAV and the vehicle are similar, and thus the birth rate and death rate for clusters around the Tx and Rx are assumed to be equal, i.e., ${\lambda }_{\mathrm{B}}=10$, ${\lambda }_{\mathrm{D}}=2$, $D_{\mathrm{c}1/\mathrm{c}2/\mathrm{c}3}=50$, $D_{\mathrm{a}1/\mathrm{a}2}=50$.

Time evolution and angular PSDs of single clusters and twin clusters are depicted in Figure 5. As mentioned before, twin clusters are randomly formed based on the simultaneous generation of new single clusters with probability $P_{twin}$. When $P_{twin}=0$, shown in (d–f), twin clusters do not have the conditions to produce. Therefore, twin clusters only contain two states, survival and death. In contrast, when $P_{twin}=0.8$ as shown in (a–c), new single clusters have a chance of merging into twin clusters, which may lead to a slower growth trend for the generation of single clusters. In this way, a suitable $P_{twin}$ can regulate the evolution speed as well as the trend of different types of clusters, and the relationship between different clusters is established. In addition, by combining with the transition region, the normalized angular PSDs are illustrated in (g–i). The birth-and-death properties are also represented in the angular domain. Moreover, the angular PSD of clusters around the Tx has more noticeable trajectory changes than clusters around the Rx due to faster angular variations resulting from closer distances. Additionally, on account of the largest delay caused by twin clusters, verified in (i), they contain the weakest power.

Figure 6 shows us the temporal ACFs for different airspeeds and at different time points. Obviously, the UAV with a high speed can generate a small temporal ACF, which can also be viewed in [40]. This is because a faster speed will cause more drastic changes in channel parameters in low-dynamic-channel environments and enhance the channel non-stationarity. In this case, a high speed will lead to a small temporal correlation. By comparing with the temporal ACFs at different time points, it can be seen that the temporal correlation changes over time, which verifies the channel non-stationarity in the time domain.

The temporal ACFs for different carrier frequencies and $k_R$ are depicted in Figure 7. Obviously, a large carrier frequency will lead to a small temporal correlation with a small coherence time, which is consistent with the conclusion in [53]. The phenomenon can be explained by the fact that a larger frequency has a smaller wavelength and leads to a more dispersive distribution of Doppler frequencies. In addition, as the concentration parameter $k_R$ increases, the temporal ACF shows a slower downtrend. $k_R$ controls the density of scatterers within clusters and scatterers are more concentrated around the average angles with the increase in $k_R$, which results in a smaller angular spread.

Figure 8 depicts the temporal ACFs for a different number of single clusters and twin clusters. It is proved that twin clusters can accelerate the rate of the curve’s decline. As is mentioned, twin clusters can be considered as combinations of single clusters. Therefore, the AoDs and AoAs are mutually independent for twin clusters and the effect leads to a more flexible spatial distribution, which will be accompanied by channel parameters with a larger distribution range, including a more dispersive Doppler spread and delay spread. Another interesting observation is that the high power of the LoS component can enhance the temporal correlation markedly. This is because the LoS component contains dominant power and effectively alleviates the channel non-stationarity in the time domain.

The temporal ACF for a single—bounced component via $C_{\mathrm{n}}^T$ and $C_{\mathrm{n}}^R$, respectively, is shown in Figure 9. The impact of the motion of the Tx, the Rx, and clusters on the temporal ACF is presented visually, and it is found that temporal ACFs dynamically change in the angular domain. This can verify the channel non-stationarity in the time domain once again. In addition, similar to [29], the temporal ACF is derived dependent on the distribution of scatterers and the area with dense scatterers results in a large temporal correlation. We also find that the horizontal movement of clusters will result in significant changes in azimuth angles. This is because the Tx, the Rx, and clusters move horizontally with different azimuth angles and thus azimuth angles contain more dramatic time-variant properties than elevation angles.

A comparison of the multi–cylinder model in [15], sphere model in [17], two–cylinder model in [19], and proposed channel model is shown in Figure 10. It can be seen that temporal ACFs share similar downward trends, which indicates the reasonability of the proposed model. In addition, we can perceive that although there are significant differences in the distribution of scatterers among models, a similar temporal correlation may be derived.

The space CCF $R_{(1,1)(1,1),(2,2)(2,2)}({\delta }_{\mathrm{T}},t)$ between the link from the $(1,1)$-th transmitting antenna to the $(1,1)$-th receiving antenna and the link from the $(2,2)$-th transmitting antenna to the $(2,2)$-th receiving antenna is explored in Figure 11. The space CCF has a faster descent speed with time, which proves the channel non-stationarity in the space domain. This is due to a greater angular spread resulting from the time-variant channel environment. Additionally, as $k_R$ increases, the space CCF shows a slower decline. This is because more dense scatterers distributed within a cluster will lead to a smaller angular spread and a smaller variance of distance differences.

Due to different spatial properties of the UPA compared to the ULA, the space CCF $R_{(1,1)(1,1),(1,{\mathrm{p}}^{\prime })(1,1)}\left({\delta }_{\mathrm{T}}\right)$ and $R_{(1,1)(1,1),({\mathrm{p}}^{\prime },1)(1,1)}\left({\delta }_{\mathrm{T}}\right)$ are explored in Figure 12. It can be seen that antennas in the same row and in the same column show different space CCFs. This can be easily explained by the fact that the UPAs can be regarded as a combination of ULAs. Thus, antennas in the same row and in the same column can be regarded as linear antennas with different space angles. As we know, the spatial angle has an influence on the space CCF due to its impact on the angular spread [19]. In addition, when two antennas are further apart, the space correlation becomes smaller. It is also influenced by a larger angular spread. Moreover, the increase in $|p-p^{\prime }|$ will result in fewer clusters shared by both antennas, which is also a factor that changes the angular spread. Therefore, the spatial correlation does not decrease smoothly and a greater spatial correlation may emerge as the antenna index increases.

Figure 13 shows time-variant characteristics of the FCF. The simultaneous effect of the continuous and random occurrence and disappearance of clusters as well as the movement of surviving clusters, the UAV, and the vehicle results in severe frequency non-stationarity. Therefore, the FCF does not change smoothly over time and the BD properties contribute to strongly random changes in the FCF.

From Figure 14, we find that the carrier frequency and airspeed have a significant influence on the Doppler PSD. As the frequency increases, the Doppler PSD becomes wider with a larger Doppler spread. In addition, the peak moves to the right. This can be explained by the fact that a small wavelength will lead to a dispersive distribution of Doppler frequencies. Similarly, the increase in the airspeed can also effectively generate a wider Doppler PSD and maintain a similar curved shape. This indicates that the velocity of the UAV will change the values of Doppler frequencies; nevertheless, it will not have a significant impact on the distribution characteristics. Additionally, the movement of the peak is also a phenomenon resulting from the enlargement of Doppler frequencies.

Figure 15 compares the RMS—DSs (a) and Doppler spreads (b) of single clusters and twin clusters, respectively. It can be found that the RMS—DS of twin clusters is the largest, which indicates that twin clusters will generate a small correlation bandwidth. This is because multi–bounced components will result in greater delays. Another observation is that the RMS—DS of twin clusters has a dispersive distribution. This implies that the emergences and movements of twin clusters are more flexible than single clusters in 3D space. Futhermore, the flexible spatial characteristics of twin clusters can also lead to a dispersive distribution of Doppler frequencies. When the AoAs or AoDs are determined by the locations of single clusters, the AoDs or AoAs of single clusters are fixed values. By contrast, the AoAs and AoDs of twin clusters are separately determined by two different single clusters randomly generated in 3D space. As a result, twin clusters have a wider angular distribution, which also enlarges the distribution range of Doppler frequencies.

6. Conclusions

A novel cluster-based low-altitude U2V non-stationary channel model with UPAs is proposed in the paper. The proposed channel model takes into account both single clusters and twin clusters. In order to model the U2V non-stationarity, an improved cluster evolution algorithm is developed to integrate continuous time evolution and array evolution. In the time evolution, new single clusters evolve based on the BD property and the formation of twin clusters relies on the simultaneous generation of new single clusters around the Tx and Rx with a correlation probability $P_{twin}$. $P_{twin}$ is proven to effectively regulate the evolutionary trends of different types of clusters. In addition, the developed array evolution algorithm for new clusters and survival clusters applied to UPAs is proved to satisfy the coherence and regionality of antenna sets. At the same time, time-variant channel parameters and channel statistical properties are derived. The simulated results prove that the proposed model can effectively capture the non-stationarity in the time, space, and frequency domain. The numerical analysis validates that twin clusters flexibly generated in 3D space contain low power and that the flexibility leads to a large distribution range of the RMS-DS and Doppler spread. Moreover, the wide distribution range of twin clusters also contributes to the fast decline in the temporal ACF as well as a large RMS-DS. We also find that the increase in the airspeed leads to a small temporal correlation and a large Doppler PSD. As scatterers become more concentrated, the space CCF and temporal ACF get large. Additionally, a large Doppler spread and temporal ACF can be observed with a large carrier frequency. The reasonability of the proposed model is also verified through comparison with other GBSMs. Undoubtedly, the speed of the UAV, the selection of carrier frequency, and changes in the scattering environment can have a perceptible effect on the channel characteristics. The characteristic difference between SB and MB components is also simulated and evaluated by single and twin clusters. These conclusions are helpful in the performance evaluation and optimization of UAV communication systems.