#### Abstract

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. This 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. The 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 time-variant Doppler power spectrum density.

#### 1. 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–7], air-to-air (A2A) channel models [8–11], and satellite-to-aircraft (S2A) channel models [12–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–10, 15–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].

##### 1.1. 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. The 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]. These 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–31]. Whereas, the distributions of path attenuation and time delay based on scatterer distribution were not analyzed in these literatures.

##### 1.2. Motivation and Contributions of This Paper

The 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. The contributions of this paper are as follows.(1)Three 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.

The 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. The superscript denotes the transpose operator. The superscript denotes the conjugate operator. is the norm of . denotes expectation operation.

#### 2. The Geometrical Structure and the Assumptions

##### 2.1. The Geometrical Structure

The geometrical structure of the presented channel model is as shown in Figure 1, where is the aircraft transmitter, is the aircraft receiver, and also the origin of the Cartesian coordinates. The ground specular reflection point is . And the diffuse scattering is assumed to happen on the flat Earth surface with the diffuse scatterer denoted with . 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, -axis is determined consequently. and are the heights of and . The time-variant coordinates of the transmitter, the receiver, the specular reflection point, and the scatterer are denoted as , , , and , respectively, and and . This 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. Thus, 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.

Based on Figure 1, the lengths of the LOS path and the specular path can be calculated by (1) and (2):

The scattering path contains subpaths, so when we discuss the subpaths, the coordinate of the th diffused scatterer is denoted as ; thus, the length of the diffuse scattering subpath can be calculated as in (3). In the following derivation, we also denote it as :

##### 2.2. The Modeling of Ground Scattering Area

The ground scatterers are assumed to be randomly distributed in the elliptic area determined by the maximum time-delay ellipsoid truncated by the XOY plane. The maximum time-delay ellipsoid is determined by (4) with and as its foci; thus, the boundary of this elliptic area is determined by (5):where 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 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 havewhere is the area of the elliptic area at the moment . The time varying of is due to the movement of the transmitter and the receiver.

##### 2.3. The Scenarios

In this paper, we define three scenarios, i.e., the GS, the VPS, and the SAS. The 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. They are as shown in Figures 2 and 3, respectively.

In the VPS, we have , , and . And in the SAS, we have .

In VPS, the scattering area is a circle, so we have , where is calculated withwhere is also assumed to be a constant for a certain system in a certain scenario.

In the GS, the area can be calculated with (8) at the moment :where and are calculated as (9) and (10), respectively:where

and are calculated as

With (8) to (15), we havewhere

In the SAS, since we have , (18) and (19) can be simplified as (20) and (21), respectively:

Thus, the PDF of the scatterers can be denoted as (22) for UD of the scatters in the abovementioned scenarios:

#### 3. 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 bywhere is the impulse response of the LOS path, is the impulse response of the specular reflection path, and 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 aswhere , with the assumption that the propagation loss attenuation factor equals to 2. is the wave length, and and are the antenna gains of the transmitter and the receiver at the moment . , where is the length of the LOS path as in (1) at the moment . in (24) is calculated aswhere , , is the velocity vector of the transmitter, , is the velocity vector of the receiver, and . The unit directional vector of the LOS path from the transmitter towards the receiver is . The impulse response of the specular reflection path is given bywhere . and can be calculated by (27) and (28):where and are the unit directional vectors towards *S* observed from and at the moment , respectively. And 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 ; 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 th subpath is calculated by (30) in accord with [21, 35]. In (30), and are the distance between the transmitter and the scatterer, and the distance between the scatterer and the receiver, respectively:

Thus, in (29), containing the phase and frequency offset in the impulse response of the diffuse scattering components, can be denoted as (31):where and are the unit directional vectors towards the th diffuse scatterer observed from and at the moment , respectively. Thus, 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):

The time-variant PSD of the scattering path can be obtained through Fourier transform with (33) as (36), where denotes :

##### 3.2. 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 , the length of the scattering path can be denoted as (37), where is the radius of the scattering area in VPS at the moment :

Thus, we have

At this moment, the Jacobian [36] can be calculated aswhere

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 for simplification of the expressions, all the parameters are time variant:

Thus, the CDF of path delay can be calculated with (42):

For the most general scenario, i.e., GS, the CDF of is calculated with (43)–(48) and can be simplified with for SAS:

Since there are no closed-form expressions of the distributions of can be obtained with the Radar function is adopted, we perform simulation and curve fitting to obtain the CDFs of 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. The distributions of the AOAs have been presented in [37].

#### 4. Simulation Assumptions and Results

##### 4.1. 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 s. At the moment 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. The 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 s, we assume that the transmitter and the receiver fly towards each other with the velocity of the transmitter m/s in Cartesian coordinates and the velocity of the receiver m/s. Thus, at the moment s, the scenario becomes VPS described in Section 3. At the moment s, we assume that m/s and m/s. Finally, at the moment s, the aircrafts become a SAS and still with the same velocities. The dielectric coefficients are set in accord with [33, 34].

##### 4.2. 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 s and 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 s and s for each moment, respectively, and the attenuation of each path are defined as normalized attenuation to the LOS path at the moment 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. This is in accord with the simulation which assumes that the scatterers are around the specular reflection point. The 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.

The time-frequency TF of the diffuse scattering path during is as shown in Figure 5, which depicts the time-variant feature of .

The normalized ACFs of the scattering path at the moments of and are shown in Figure 6.

From Figure 6, we can observe that the main lobe of the normalized ACF at is much wider than that at . This depicts different correlation time periods. We can also observe that the fluctuations of the former curve are not so obvious as the latter. These show the nonstationary feature of .

The normalized Doppler PSDs at the moments s and 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 s, when the transmitter and the receivers are performing a vertical pass-by.

##### 4.3. Simulations of the Distributions of Path Delay, Path Attenuation, and Normalized Doppler Shift of the Scattering Path

The CDFs of at the moments s, s, and s are, as shown in Figure 8, with .

From Figure 8, we can see that CDF of is apparently time variant. The simulated CDFs coincide with the theoretical results well. With similar geometrical structures, similar features can be observed between the curves at the moments of s and s.

The CDFs of at the moments s, s, and s are as shown in Figure 9 with .

From Figure 9, the nonstationary feature of the channel can be clearly observed, as well as the similarity between the CDFs at s and s with similar geometrical structures.

For better observation of the distribution of , we present the PDF of in GS in Figure 10 and fit it with the 7th order Gaussian function, as (49), where , , and . None of the well-known distributions such as Gaussian distribution, Rice distribution, and Rayleigh distribution can fit this distribution of well. And the 7th order Gaussian function is much easier to calculate than the well-known Suzuki distribution [39].

The PDF of in VPS is illustrated in Figure 11, and fitted with two terms’ exponential function as (50). Compared with the GS and SAS, the simpler geometrical structure of the VPS leads to simpler expression of the PDF of .

The simulated PDF of in SAS, i.e., at the moment s, is presented in Figure 12 and is fitted with the 4th order Gaussian function, as (51), where , . Compared with the GS, a simpler expression of the PDF of is obtained with best goodness of fit.

#### 5. 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. Then, 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.

The 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–42].

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.