Analytical RS-GBSM-Based Nonstationary Low-Altitude Air-to-Air Channel Modeling over Open Area

Aerial communication is very flexible due to almost no restrictions on geographical conditions. In recent years, with the development and application of the unmanned aerial vehicle, the air-to-air communication attracts dense interests from the researchers. More accurate and precise channel modeling for air-to-air communication is a new hot topic because of its essential role in the performance evaluation of the systems.*is paper presents an analytical nonstationary regular-shaped geometry-based statistical model for low-altitude air-to-air communication over an open area with considerations on ground scattering. Analytical expressions of the channel impulse response and the autocorrelation functions based on the three-ray model are derived. Based on the assumption of uniform distribution of the ground scatterers, the distributions of the channel coefficients such as time delay and path attenuation are derived, simulated, compared, and fitted. *e nonstationary characteristics of the channel are observed through the time-variant distributions of the channel coefficients as well as the time-variant autocorrelated functions and timevariant Doppler power spectrum density.


Introduction
Aerial communication is useful either as a compensation of the land mobile communication network or as an independent mobile network. It is very flexible and convenient due to almost no restriction on geographical conditions. As an essential part of the performance evaluation of the aerial communication systems, channel modeling of the aerial communication system is a very hot topic [1,2]. In the past years, aeronautical channel models, including air-to-ground (A2G) channel models [3][4][5][6][7], airto-air (A2A) channel models [8][9][10][11], and satellite-toaircraft (S2A) channel models [12][13][14], have been widely discussed and developed. As emerging applications of unmanned aerial vehicle (UAV) communications, precise low-altitude A2A channel modeling attracts interests from the researchers, as in [8][9][10][15][16][17][18]; measurements were also performed in [19,20]. However, there are still lots of work to do in this area, for example, to obtain the analytical expressions of the nonstationary channel with different geometrical structures and to describe various channel characteristics due to the distributions of the scatterers [16,17].

Related Works.
In 1973, Bello [21] proposed the methodologies of modeling satellite-to-aircraft channel; after that, the three-path model proposed in [21], i.e., the line of sight (LOS) path, the specular reflection path, and the diffuse scattering path, has been widely adopted in the modeling of mobile-to-mobile (M2M) communication channels including A2G channel, vehicle-to-vehicle (V2V) channel, and A2A channel; the specular reflection path was combined with the scattering paths in recent literature. e way of modeling the M2M channel models, i.e., formulating the channel impulse response (CIR) of different paths with terms of time delay, path attenuation, and Doppler shift, has also been widely adopted. In recent years, with more precise considerations on the movements in the channel, nonstationary characteristics have been widely discussed in V2V channel modeling [22,23] and also have been developed for A2G scenario [6,7] and A2A scenario [9,16,17]. ese models mostly belonged to geometry-based statistical models (GBSM).
On the modeling of the geometrical structures, Lian et al. [24] used two-dimensional (2D) one-ring and three-dimensional (3D) multicylinders to model moving scatterers and near and far stationary scatterers for high-altitude platform. Zhang and Cheng [25] proposed a two-cylinder model. In [6], as a typical example of regular-shaped GBSM (RS-GBSM), a 3D ellipsoidal model was proposed for A2G communication, and Zhang et al. [17] expanded the 3D ellipsoidal model to A2A communication.
On the assumed distributions of the scatterers and subsequent derivations, Gulfam et al. [6] derived the distributions of the received power and Doppler shift and simulated the distributions of the azimuth and elevation angles of arrival (AOAs) based on the uniform distribution (UD) of the scatterers; Zhu et al. [26] adopted a 3D von Mises Fisher (VMF) AOA distribution scenario. Zhang et al. [17] mentioned the distribution of received power caused by the stochastical characteristics of the scatterers, but the authors have not presented a more precise path attenuation calculation.
Besides, measurements were performed in [19,20]. Liu et al. [20] presented a measurement-based A2A channel model with an illustration of the distribution of shadow fading, but further discussions on other channel coefficients are expected. In [9,27,28], ground scattering was modeled and the Doppler effect was simulated in detail. In [28], a simpler closed-form expression of Doppler shift was obtained through utilizing the prolate spheroidal coordinate system, and further development could be found in [29][30][31]. Whereas, the distributions of path attenuation and time delay based on scatterer distribution were not analyzed in these literatures.

Motivation and Contributions of is Paper.
e motivation of this paper is to develop an analytical nonstationary RS-GBSM for low-altitude A2A communications over an open area, with the scatterers assumed to obey UD in a maximum delay-determined ellipsoidal scattering region. In this paper, the formulation of the channel and the derivation of the time-variant distributions of the channel coefficients are presented. e contributions of this paper are as follows.
(1) ree paths, i.e., the LOS path, the ground specular reflection path, and the ground scattering path, are modeled in this paper. Analytical expressions of CIR, transfer function (TF), autocorrelation function (ACF), and Doppler power spectrum density (PSD), are derived. (2) Nonstationary characteristics of the channel are observed through the CIRs, TFs, ACFs, Doppler PSD, etc.
(3) Closed-form expressions of the distributions of time delay are derived in three scenarios, i.e., the general scenario (GS), the vertical pass-by scenario (VPS), and the same altitude scenario (SAS). (4) Distributions of path attenuation based on radar function in these three scenarios are simulated and fitted.
e rest of this paper is organized as follows. In Section 2, we illustrate the geometrical structure and assumptions. In Section 3, we develop the nonstationary A2A channel model, formulate the analytical CIR and the time-variant TFs, and derive the time-frequency ACFs. In Section 4, we deduce the cumulative distribution functions (CDFs) of time delay and discuss the distributions of path attenuation in the GS, the VPS, and the SAS. In Section 4.2, we present the simulation results, fitting functions, and analyze them. Finally, in Section 5, we provide conclusions.
Notation: boldface uppercase letters stand for matrices while boldface lowercase letters represent vectors. e superscript (·) T denotes the transpose operator. e superscript (·) * denotes the conjugate operator. |x| is the norm of x. E · { } denotes expectation operation.

e Geometrical Structure.
e geometrical structure of the presented channel model is as shown in Figure 1, where T is the aircraft transmitter, R is the aircraft receiver, and also the origin O of the Cartesian coordinates. e ground specular reflection point is S. And the diffuse scattering is assumed to happen on the flat Earth surface with the diffuse scatterer denoted with D. D is distributed in a maximum delay-determined area which is determined in Section 2.2. Major parameters used in this paper are listed in Table 1.
Triangle RST determines the YOZ plane; thus, x-axis is determined consequently. h T and h R are the heights of T and R. e time-variant coordinates of the transmitter, the receiver, the specular reflection point, and the scatterer are denoted as . is general geometrical structure illustrated in Figure 1 can be expanded to the scenarios in [6,17] conveniently. In Figure 1, since the transmitter and the receiver are moving during communication, the distances and the AOAs are time variant. us, in most of the situations, the channel is highly nonstationary. According to the general method of developing aeronautical channel models, such as presented in [5,21], three paths, including the line-of-sight path, the reflection path, and the diffuse scattering paths, are considered in this paper and are assumed to be uncorrelated to each other.

2
International Journal of Antennas and Propagation Based on Figure 1, the lengths of the LOS path and the specular path can be calculated by (1) and (2): (2) e scattering path contains subpaths, so when we discuss the subpaths, the coordinate of the lth diffused scatterer is denoted as (x D,l , y D,l , 0); thus, the length of the diffuse scattering subpath can be calculated as in (3). In the following derivation, we also denote it as r DIF :

e Modeling of Ground Scattering
Area. e ground scatterers are assumed to be randomly distributed in the elliptic area determined by the maximum time-delay ellipsoid truncated by the XOY plane. e maximum timedelay ellipsoid is determined by (4) with T and R as its foci; thus, the boundary of this elliptic area is determined by (5): where τ max is the maximum delay of the scattering paths in second and c is the speed of light. Compared with the delay of the LOS path, we define τ max � ϵτ LOS as the maximum delay of the diffuse scattering path that can be detected, which assumed to be a constant when a certain communication system is under consideration. Since we assume that the scatterers obey UD in the elliptic area, we have International Journal of Antennas and Propagation e time-variant 3D coordinate of the lth scatterer Φ D e azimuth angle between the scatterer D and y-axis in VPS α LOS (t) Path attenuation of the LOS path α SPE (t) Path attenuation of the specular reflection path α DIF (t) Path attenuation of the scattering path α DIF,l (t) Path attenuation of the lth scattering path Δt Time interval Internal variables 4 International Journal of Antennas and Propagation where S E (t) is the area of the elliptic area at the moment t. e time varying of S E (t) is due to the movement of the transmitter and the receiver.

e Scenarios.
In this paper, we define three scenarios, i.e., the GS, the VPS, and the SAS. e GS is illustrated in Figure 1, the VPS stands for the scenario that the transmitter and the receiver have the same XOY coordinates but different altitudes, and the SAS stands for the scenario that the transmitter and the receiver are with the same altitude but different XOY coordinates. ey are as shown in Figures 2  and 3, respectively.
In the VPS, we have In VPS, the scattering area is a circle, so we have where r DIF max � cτ max is also assumed to be a constant for a certain system in a certain scenario.
In the GS, the area S E can be calculated with (8) at the moment t: where y D min (t) and y D max (t) are calculated as (9) and (10), respectively: where x D upper and x D lower are calculated as x International Journal of Antennas and Propagation 5 With (8) to (15), we have In the SAS, since we have h T � h R , (18) and (19) can be simplified as (20) and (21), respectively: us, the PDF of the scatterers can be denoted as (22) for UD of the scatters in the abovementioned scenarios:

The Analytical Nonstationary A2A Channel Model
3.1. Channel Formulation. As mentioned above, in accord with [5,21], three-path model is adopted in this paper; thus, the CIR of the presented channel model can be given by where h LOS (t, τ) is the impulse response of the LOS path, h SPE (t, τ) is the impulse response of the specular reflection path, and h DIF (t, τ) is the impulse response of the diffuse scattering path. By considering the large-scale propagation attenuation, propagation delay, and Doppler shift, the impulse response of the LOS path is denoted as where α LOS (t) � λ( , with the assumption that the propagation loss attenuation factor equals to 2. λ is the wave length, and G LOS T (t) and G LOS R (t) are the antenna gains of the transmitter and the receiver at the moment t. τ LOS (t) � (d LOS (t)/c), where d LOS (t) is the length of the LOS path as in (1) at the moment t. g LOS (t) in (24) is calculated as International Journal of Antennas and Propagation e unit directional vector of the LOS path from the transmitter towards the receiver is n RT (t) � (n RT x(t), n RT y(t), n RT z(t)) T . e impulse response of the specular reflection path is given by where τ SPE � d SPE (t)/c. α SPE (t) and g SPE (t) can be calculated by (27) and (28): where n TS and n RS are the unit directional vectors towards S observed from T and R at the moment t, respectively. And Γ(t) is the reflection coefficient where the dielectric property of the reflection surface is calculated as in [32,33] and [34].
In the simulation of this paper, we assume that the scattering surface is the sea surface. Suppose that the number of the scatterers is L; then, the impulse response of the scattering path is as When the radar cross section (RCS) σ is constant for each scatterer, the amplitude attenuation of the lth subpath α l (t) is calculated by (30) in accord with [21,35]. In (30), d TD,l (t) and d DR,l (t) are the distance between the transmitter and the scatterer, and the distance between the scatterer and the receiver, respectively: us, g DIF,l (t) in (29), containing the phase and frequency offset in the impulse response of the diffuse scattering components, can be denoted as (31): where n TD,l (t) � (n TD,l x (t), n TD,l y (t), n TD,l z (t)) T and n RD,l (t) � (n RD,l x (t), n RD,l y (t), n RD,l z (t)) T are the unit directional vectors towards the lth diffuse scatterer D observed from T and R at the moment t, respectively. us, the TF and ACF of the scattering paths are calculated as (32) and (33) With uncorrelated assumption, the TF and ACF of the scattering paths can be stationary in frequency domain; thus, (32) and (33) become (34) and (35): e time-variant PSD of the scattering path can be obtained through Fourier transform with (33) as (36), where u denotes Δt:

Distributions of Channel Coefficients of Scattering Path.
In VPS, the distribution of the scattering path delay can be derived as follows.
As shown in Figure 2, at the moment t, the length of the scattering path r DIF (t) can be denoted as (37), where r D (t) is the radius of the scattering area in VPS at the moment t: us, we have . (38) At this moment, the Jacobian [36] can be calculated as International Journal of Antennas and Propagation where ξ � According to the UD of the scatterers, we also can denote the PDF of the scatterers as (41). It is to be mentioned that, although we omit the variable t for simplification of the expressions, all the parameters are time variant: us, the CDF of path delay τ DIF can be calculated with (42):

International Journal of Antennas and Propagation
For the most general scenario, i.e., GS, the CDF of τ DIF is calculated with (43)-(48) and can be simplified with h T � h R for SAS: Since there are no closed-form expressions of the distributions of α DIF can be obtained with the Radar function is adopted, we perform simulation and curve fitting to obtain the CDFs of α DIF in GS, VPS, and SAS. Furthermore, since the PDF of the Doppler shift is proportional to the Doppler PSD and also has been derived in [28]; we do not derive the distribution of Doppler shift in this paper. e distributions of the AOAs have been presented in [37].

Simulation Assumptions.
According to similar scenarios in [5,9] and [11], the simulation parameters are set and listed in Table 2. In our simulation, we assume that the simulation starts at the moment t � 0 s. At the moment t � 0 s, the transmitter and the receiver are at different altitudes, which are corresponding to the GS, e.g., the altitude of the transmitter is 305 m and the altitude of the receiver is 610 m. e initial distance between the transmitter and the receiver is assumed to be 745 m, i.e., the ground distance is 680 m. At the moment t � 0 s, we assume that the transmitter and the receiver fly towards each other with the velocity of the transmitter v T � (0, 68, 0) T m/s in Cartesian coordinates and the velocity of the receiver v R � (0, −68, 0) T m/s. us, at the moment t � 5 s, the scenario becomes VPS described in Section 3. At the moment t � 5 s, we assume that v T � (0, −68, 0) T m/s and v R � (0, 68, −61) T m/s. Finally, at the moment t � 10 s, the aircrafts become a SAS and still with the same velocities. e dielectric coefficients are set in accord with [33,34].

Simulation
Results of CIR, TF, ACF, and PSD. In our simulations of the CIRs, the tapped delay lines model is adopted. Figure 4 shows the impulse responses of the presented model at the moments t � 0 s and t � 5 s. Fifty taps are illustrated in this figure with equal delay paths overlapped. For convenience, the tap delays in Figure 4 are set as the time delay to t � 0 s and t � 5 s for each moment, respectively, and the attenuation of each path are defined as normalized attenuation to the LOS path at the moment t � 0 s. It can be seen from Figure 4 that the LOS path, which is also the path with minimum delay, can be distinguished from other paths easily, while the specular reflection path is adjacent to the diffuse scatter paths. is is in accord with the simulation which assumes that the scatterers are around the specular reflection point. e time-varying feature of the CIRs is very obvious in Figure 4. It is to be noted that, in Figure 4, each diffuse component is much smaller than the LOS component, but it does not mean that we should neglect the diffuse component. Modeling the diffuse components is the most important part of channel modeling since they introduce stochastics of the channel. For the realization of computer simulation, we take a finite number of the diffuse scattering paths. However, in many pieces of [11,18,38], the number of diffuse paths goes to infinity. With this assumption of infinite diffuse paths, the amplitude of each scattering tap (distinguishable diffuse paths) will become much higher. e time-frequency TF of the diffuse scattering path during t ∈ (0, 0.12)s is as shown in Figure 5, which depicts the time-variant feature of H DIF (t, f). e normalized ACFs of the scattering path at the moments of t � 0s and t � 5s are shown in Figure 6.
From Figure 6, we can observe that the main lobe of the normalized ACF at t � 0s is much wider than that at t � 5s.
is depicts different correlation time periods. We can also observe that the fluctuations of the former curve are not so obvious as the latter. ese show the nonstationary feature of H DIF (t, f). e normalized Doppler PSDs at the moments t � 0 s and t � 5 s are shown in Figure 7. Figure 7 shows the decreasing of the maximum normalized Doppler in accord with the movements of the transmitter and the receiver towards each other with the increasing of elevation angles. Furthermore, the Doppler PSD is a bilateral spectrum because at the moment around t � 5 s, when the transmitter and the receivers are performing a vertical pass-by.

Simulations of the Distributions of Path Delay, Path Attenuation, and Normalized Doppler Shift of the Scattering
Path.
e CDFs of τ DIF at the moments t � 0 s, t � 5 s, and t � 10 s are, as shown in Figure 8, with ϵ � 3.57.
From Figure 8, we can see that CDF of τ DIF is apparently time variant. e simulated CDFs coincide with the theoretical results well. With similar geometrical structures, similar features can be observed between the curves at the moments of t � 0 s and t � 10 s. e CDFs of α DIF at the moments t � 0 s, t � 5 s, and t � 10 s are as shown in Figure 9 with ϵ � 3.57.
From Figure 9, the nonstationary feature of the channel can be clearly observed, as well as the similarity between the CDFs at t � 0 s and t � 10 s with similar geometrical structures.
For better observation of the distribution of α DIF , we present the PDF of α DIF in GS in Figure 10 and fit it with the 7th order Gaussian function, as (49), where a 1 � 1.858E − 07, and a 7 � 1.581E + 05, b 7 � 1.362E − 06, c7 � 2.063E − 06. None of the well-known distributions such as Gaussian distribution, Rice distribution, and Rayleigh distribution can fit this distribution of α DIF well. And the 7th order Gaussian function is much easier to calculate than the well-known Suzuki distribution [39].

Conclusion
We present an analytical RS-GBSM-based channel model for low-altitude A2A communication over an open area in this paper. Based on the assumption of the scatterers distributed in a maximum-delay-determined scattering area, we formulate the CIRs for the LOS path, the specular-reflected path, and the scattering path. Furthermore, we derived the TF, the ACF, and the Doppler PSD of the scattering path since the scattering path is the most complicated path. en, with the UD of the scatterers, we derive the closed-form expression of the distributions of the path delay of the scattering path in VPS, SAS, and GS. Simulation results show a good coincidence between the theoretical CDFs and the simulated CDFs. Since the closed-form expression of the path attenuation distributions based on the radar function is very hard to be derived, simulation and curve fitting are adopted, and we fit these functions and present the corresponding parameters. e work performed in this paper is one of the trials towards the analytical nonstationary channel modeling based on the geometrical distribution of the scatterers. In further research, we will consider more kinds of distributions. For more comprehensive structures, i.e., including the modeling of local scatterers around the ground station or the aircraft, if necessary, we can refer to [40][41][42].

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.