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

Authors: Zaixue Wei; Qipeng Tang
Journal: International Journal of Antennas and Propagation
 (2022)

---

## 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.


---

## Body

## 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 Section2, 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⋅T denotes the transpose operator. The superscript ⋅∗ denotes the conjugate operator. x is the norm of x. E⋅ denotes expectation operation.


## 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 Section2, 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⋅T denotes the transpose operator. The superscript ⋅∗ denotes the conjugate operator. x is the norm of x. E⋅ 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 Figure1, where T is the aircraft transmitter, R is the aircraft receiver, and also the origin O of the Cartesian coordinates. The 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.Figure 1
The geometrical structure of the system model with flat Earth surface as the scattering surface.Table 1
Major parameters and symbols this paper.
SymbolsDefinitioncThe speed of lightDThe scattererdLOStTime variant distance between the transmitter and the receiver (length of the LOS path)dSPEtThe time-variant length of the specular reflection pathdTDtThe time-variant distance between the transmitterT and the scatterer DdRDtThe time-variant distance between the specular reflection pointD and the receiver RdDIF,ltThe time-variant length of thelth scattering subpathfFrequencyfcCentral frequencyGTLOStThe antenna gain of the receiver at the momentt along the LOS pathGRLOStThe antenna gain of the receiver at the momentt along the LOS pathGTSPEtThe antenna gain of the receiver at the momentt along the specular reflection pathGRSPEtThe antenna gain of the receiver at the momentt along the specular reflection pathGT,lDIFtThe antenna gain of the receiver at the momentt along the specular reflection pathGR,lDIFtThe antenna gain of the receiver at the momentt along the specular reflection pathhRThe altitude of the receiverhRtThe time-variant altitude of the receiverhTThe altitude of the transmitterhTtThe time-variant altitude of the transmitterht,τCIR of the channelhLOSt,τCIR of the LOS pathhSPEt,τCIR of the specular reflection pathhDIFt,τCIR of the scattering pathLNumber of the scatterersnRTtThe time-variant unity directional vector of the LOS path from the transmitter towards the receivernTOtThe time-variant unit directional vector towards O observed from TnROtThe time-variant unit directional vector towards O observed from RnTD,ltThe time-variant unit directional vector towards thelth scatterer observed from the transmitter TnRD,ltThe time-variant unit directional vector towards thel th scatterer observed from the receiver ROThe origin of the XOY plane, also the specular reflection pointrDThe radius of the scattering area in VPSrDtThe time-variant radius of the scattering area in VPSrDIFThe length of the scattering pathrDIFtThe time-variant length of the scattering pathRThe receiverSThe specular reflection pointSEtThe time-variant area of the scattering ellipsetTimeTThe transmittervRtTime variant velocity vector of the receivervTtTime variant velocity vector of the transmitterxTt,yTt,zTtThe time-variant 3D coordinate of the transmitterxRt,yRt,zRtThe time-variant 3D coordinate of the receiverxDt,yDt,0tThe time-variant 3D coordinate of the scattererxD,lt,yD,lt,0tThe time-variant 3D coordinate of thelth scattererΦDThe azimuth angle between the scattererD and y-axis in VPSαLOStPath attenuation of the LOS pathαSPEtPath attenuation of the specular reflection pathαDIFtPath attenuation of the scattering pathαDIF,ltPath attenuation of thelth scattering pathΔtTime intervalΔfFrequency intervalΓtTime variant reflection coefficientγDIFNormalized Doppler shift of the scattering pathσRadar cross section (RCS)τDelayτLOStTime variant path delay of the LOS path at the momenttτSPEtTime variant path delay of the specular reflection path at the momenttτDIFPath delay of the scattering pathτDIFtTime variant path delay of the scattering path at the momenttτDIF,ltTime variant path delay of thelth scattering path at the moment tθThe reflection angleλThe wavelengthξ,ξit,i=1,…,3Internal variablesμit,i=1,…,3Internal variablesTriangle RST determines the YOZ plane; thus,x-axis is determined consequently. hT and hR are the heights of T and R. The time-variant coordinates of the transmitter, the receiver, the specular reflection point, and the scatterer are denoted as TxTt,yTt,zTt, R0,0,0, SxSt,ySt,zSt, and DxDt,yDt,zDt, respectively, and zTt=hTt−hRt and zSt=zDt=−hRt. 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 Figure1, the lengths of the LOS path and the specular path can be calculated by (1) and (2):(1)dLOSt=xT2t+yT2t+zT2t,(2)dSPEt=dTSt+dSRt=hT2t+xTt−xSt2+yTt−ySt2+hR2t+xS2t+yS2t.The scattering path contains subpaths, so when we discuss the subpaths, the coordinate of thelth diffused scatterer is denoted as xD,l,yD,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 rDIF:(3)dDIF,lt=dRD,lt+dTD,lt=xD,l2t+yD,l2t+hR2t+hT2t+xTt−xD,lt2+yTt−yD,lt2.


### 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 T and R as its foci; thus, the boundary of this elliptic area is determined by (5):(4)xD2t+yD2t+zD2t+xDt−xTt2+yDt−yTt2+zDt−zTt2=cτmax,(5)xD2t+yD2t+hR2t+xDt−xTt2+yDt−yTt2+hT2t=cτmax,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(6)fXDYDxD,yD,t=1SEt,where SEt is the area of the elliptic area at the moment t. The time varying of SEt 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 Figure1, 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.Figure 2
The geometrical structure of the VPS.Figure 3
The geometrical structure of the SAS.In the VPS, we havexTt=xRt=0, yTt=yRt=0, and zTt=hTt−hRt. And in the SAS, we have zTt=zRt=hRt.In VPS, the scattering area is a circle, so we haveSEt=1/πrDmax2t, where rDmaxt is calculated with(7)rDmaxt=12rDIFmaxhRt+hTt−rDIFmax⋅hTt−hRt−rDIFmax⋅hTt−hRt+rDIFmax⋅hRt+hTt+rDIFmax,where rDIFmax=cτmax is also assumed to be a constant for a certain system in a certain scenario.In the GS, the areaSE can be calculated with (8) at the moment t:(8)SEt=∫yDmintyDmaxt∫xD_lowertxD_uppertdxDdyD=∫yDmintyDmaxtxD_uppert−xD_lowertdyD,where yDmint and yDmaxt are calculated as (9) and (10), respectively:(9)yDmint=ξ3t−ξ1tξ2t,(10)yDmaxt=ξ3t+ξ1tξ2t,where(11)ξ1t=rDIFmax−rDIFmax2+hR2t−2hRthTt+hT2t+yT2t⋅−rDIFmax2+hR2t+2hRthTt+hT2t+yT2t,(12)ξ2t=2yT2t−2rDIFmax2,(13)ξ3t=hT2tyTt−hR2tyTt+yT3t−rDIFmax2yTt.xD_upper and xD_lower are calculated as(14)xD_uppert=12rDIFmaxrDIFmax4−2hR2trDIFmax2−2hT2trDIFmax2−4rDIFmax2yD2+4rDIFmax2yDyTt−2rDIFmax2yT2t+hR4t−2hR2thT2t+4hR2tyDyTt−2hR2tyT2t+hT4t−4hT2tyDyTt+2hT2tyT2t+4yD2yT2t−4yDyT3t+yT4t,(15)xD_lowert=−12rDIFmaxrDIFmax4−2hR2trDIFmax2−2hT2trDIFmax2−4rDIFmax2yD2+4rDIFmax2yDyTt−2rDIFmax2yT2t+hR4t−2hR2thT2t+4hR2tyDyTt−2hR2tyT2t+hT4t−4hT2tyDyTt+2hT2tyT2t+4yD2yT2t−4yDyT3t+yT4t.With (8) to (15), we have(16)SEt=1rDIFmax4μ1tyDmaxt+μ2tμ1tyDmax2t+μ2tyDmaxt+μ3t+logμ1tyDmax2t+μ2tyDmaxt+μ3t+μ2t+2μ1tyDmaxt2μ1t4μ1tμ3t−μ22t8μ13/2t−4μ1tyDmint+μ2tμ1tyDmin2t+μ2tyDmint+μ3t−logμ1tyDmin2t+μ2tyDmint+μ3t+μ2t+2μ1tyDmint2μ1t4μ1tμ3t−μ22t8μ13/2t,where(17)μ1t=−4rDIFmax2+4yT2t,(18)μ2t=4rDIFmax2yTt+4hR2tyTt−4hT2tyT−4yT3t,(19)μ3t=rDIFmax4−2hR2trDIFmax2−2hT2trDIFmax2−2rDIFmax2yT2t+hR4t−2hR2thT2t−2hR2tyT2t+hT4t+2hT2tyT2t+yT4t.In the SAS, since we havehT=hR, (18) and (19) can be simplified as (20) and (21), respectively:(20)μ2t=4rDIFmax2yTt−4yT3t,(21)μ3t=rDIFmax4−4hR2trDIFmax2−2rDIFmax2yT2t+2hR4t−2hR2thT2+yT4t.Thus, the PDF of the scatterers can be denoted as (22) for UD of the scatters in the abovementioned scenarios:(22)fXDYDxD,yD,t=1πrDmax2t,VPS,rDIFmax4μ1tyDmaxt+μ2t⋅μ1tyDmax2t+μ2tyDmaxt+μ3t+logμ1tyDmax2t+μ2tyDmaxt+μ3t+μ2t+2μ1tyDmaxt/2μ1t8μ13/2t⋅4μ1tμ3t−μ22t−4μ1tyDmint+μ2t⋅μ1tyDmin2t+μ2tyDmint+μ3t−logμ1tyDmin2t+μ2tyDmint+μ3t+μ2t+2μ1tyDmint/2μ1t8μ13/2t⋅4μ1tμ3t−μ22t−1,GS&SAS.


## 2.1. The Geometrical Structure

The geometrical structure of the presented channel model is as shown in Figure1, where T is the aircraft transmitter, R is the aircraft receiver, and also the origin O of the Cartesian coordinates. The 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.Figure 1
The geometrical structure of the system model with flat Earth surface as the scattering surface.Table 1
Major parameters and symbols this paper.
SymbolsDefinitioncThe speed of lightDThe scattererdLOStTime variant distance between the transmitter and the receiver (length of the LOS path)dSPEtThe time-variant length of the specular reflection pathdTDtThe time-variant distance between the transmitterT and the scatterer DdRDtThe time-variant distance between the specular reflection pointD and the receiver RdDIF,ltThe time-variant length of thelth scattering subpathfFrequencyfcCentral frequencyGTLOStThe antenna gain of the receiver at the momentt along the LOS pathGRLOStThe antenna gain of the receiver at the momentt along the LOS pathGTSPEtThe antenna gain of the receiver at the momentt along the specular reflection pathGRSPEtThe antenna gain of the receiver at the momentt along the specular reflection pathGT,lDIFtThe antenna gain of the receiver at the momentt along the specular reflection pathGR,lDIFtThe antenna gain of the receiver at the momentt along the specular reflection pathhRThe altitude of the receiverhRtThe time-variant altitude of the receiverhTThe altitude of the transmitterhTtThe time-variant altitude of the transmitterht,τCIR of the channelhLOSt,τCIR of the LOS pathhSPEt,τCIR of the specular reflection pathhDIFt,τCIR of the scattering pathLNumber of the scatterersnRTtThe time-variant unity directional vector of the LOS path from the transmitter towards the receivernTOtThe time-variant unit directional vector towards O observed from TnROtThe time-variant unit directional vector towards O observed from RnTD,ltThe time-variant unit directional vector towards thelth scatterer observed from the transmitter TnRD,ltThe time-variant unit directional vector towards thel th scatterer observed from the receiver ROThe origin of the XOY plane, also the specular reflection pointrDThe radius of the scattering area in VPSrDtThe time-variant radius of the scattering area in VPSrDIFThe length of the scattering pathrDIFtThe time-variant length of the scattering pathRThe receiverSThe specular reflection pointSEtThe time-variant area of the scattering ellipsetTimeTThe transmittervRtTime variant velocity vector of the receivervTtTime variant velocity vector of the transmitterxTt,yTt,zTtThe time-variant 3D coordinate of the transmitterxRt,yRt,zRtThe time-variant 3D coordinate of the receiverxDt,yDt,0tThe time-variant 3D coordinate of the scattererxD,lt,yD,lt,0tThe time-variant 3D coordinate of thelth scattererΦDThe azimuth angle between the scattererD and y-axis in VPSαLOStPath attenuation of the LOS pathαSPEtPath attenuation of the specular reflection pathαDIFtPath attenuation of the scattering pathαDIF,ltPath attenuation of thelth scattering pathΔtTime intervalΔfFrequency intervalΓtTime variant reflection coefficientγDIFNormalized Doppler shift of the scattering pathσRadar cross section (RCS)τDelayτLOStTime variant path delay of the LOS path at the momenttτSPEtTime variant path delay of the specular reflection path at the momenttτDIFPath delay of the scattering pathτDIFtTime variant path delay of the scattering path at the momenttτDIF,ltTime variant path delay of thelth scattering path at the moment tθThe reflection angleλThe wavelengthξ,ξit,i=1,…,3Internal variablesμit,i=1,…,3Internal variablesTriangle RST determines the YOZ plane; thus,x-axis is determined consequently. hT and hR are the heights of T and R. The time-variant coordinates of the transmitter, the receiver, the specular reflection point, and the scatterer are denoted as TxTt,yTt,zTt, R0,0,0, SxSt,ySt,zSt, and DxDt,yDt,zDt, respectively, and zTt=hTt−hRt and zSt=zDt=−hRt. 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 Figure1, the lengths of the LOS path and the specular path can be calculated by (1) and (2):(1)dLOSt=xT2t+yT2t+zT2t,(2)dSPEt=dTSt+dSRt=hT2t+xTt−xSt2+yTt−ySt2+hR2t+xS2t+yS2t.The scattering path contains subpaths, so when we discuss the subpaths, the coordinate of thelth diffused scatterer is denoted as xD,l,yD,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 rDIF:(3)dDIF,lt=dRD,lt+dTD,lt=xD,l2t+yD,l2t+hR2t+hT2t+xTt−xD,lt2+yTt−yD,lt2.


## 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 T and R as its foci; thus, the boundary of this elliptic area is determined by (5):(4)xD2t+yD2t+zD2t+xDt−xTt2+yDt−yTt2+zDt−zTt2=cτmax,(5)xD2t+yD2t+hR2t+xDt−xTt2+yDt−yTt2+hT2t=cτmax,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(6)fXDYDxD,yD,t=1SEt,where SEt is the area of the elliptic area at the moment t. The time varying of SEt 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 Figure1, 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.Figure 2
The geometrical structure of the VPS.Figure 3
The geometrical structure of the SAS.In the VPS, we havexTt=xRt=0, yTt=yRt=0, and zTt=hTt−hRt. And in the SAS, we have zTt=zRt=hRt.In VPS, the scattering area is a circle, so we haveSEt=1/πrDmax2t, where rDmaxt is calculated with(7)rDmaxt=12rDIFmaxhRt+hTt−rDIFmax⋅hTt−hRt−rDIFmax⋅hTt−hRt+rDIFmax⋅hRt+hTt+rDIFmax,where rDIFmax=cτmax is also assumed to be a constant for a certain system in a certain scenario.In the GS, the areaSE can be calculated with (8) at the moment t:(8)SEt=∫yDmintyDmaxt∫xD_lowertxD_uppertdxDdyD=∫yDmintyDmaxtxD_uppert−xD_lowertdyD,where yDmint and yDmaxt are calculated as (9) and (10), respectively:(9)yDmint=ξ3t−ξ1tξ2t,(10)yDmaxt=ξ3t+ξ1tξ2t,where(11)ξ1t=rDIFmax−rDIFmax2+hR2t−2hRthTt+hT2t+yT2t⋅−rDIFmax2+hR2t+2hRthTt+hT2t+yT2t,(12)ξ2t=2yT2t−2rDIFmax2,(13)ξ3t=hT2tyTt−hR2tyTt+yT3t−rDIFmax2yTt.xD_upper and xD_lower are calculated as(14)xD_uppert=12rDIFmaxrDIFmax4−2hR2trDIFmax2−2hT2trDIFmax2−4rDIFmax2yD2+4rDIFmax2yDyTt−2rDIFmax2yT2t+hR4t−2hR2thT2t+4hR2tyDyTt−2hR2tyT2t+hT4t−4hT2tyDyTt+2hT2tyT2t+4yD2yT2t−4yDyT3t+yT4t,(15)xD_lowert=−12rDIFmaxrDIFmax4−2hR2trDIFmax2−2hT2trDIFmax2−4rDIFmax2yD2+4rDIFmax2yDyTt−2rDIFmax2yT2t+hR4t−2hR2thT2t+4hR2tyDyTt−2hR2tyT2t+hT4t−4hT2tyDyTt+2hT2tyT2t+4yD2yT2t−4yDyT3t+yT4t.With (8) to (15), we have(16)SEt=1rDIFmax4μ1tyDmaxt+μ2tμ1tyDmax2t+μ2tyDmaxt+μ3t+logμ1tyDmax2t+μ2tyDmaxt+μ3t+μ2t+2μ1tyDmaxt2μ1t4μ1tμ3t−μ22t8μ13/2t−4μ1tyDmint+μ2tμ1tyDmin2t+μ2tyDmint+μ3t−logμ1tyDmin2t+μ2tyDmint+μ3t+μ2t+2μ1tyDmint2μ1t4μ1tμ3t−μ22t8μ13/2t,where(17)μ1t=−4rDIFmax2+4yT2t,(18)μ2t=4rDIFmax2yTt+4hR2tyTt−4hT2tyT−4yT3t,(19)μ3t=rDIFmax4−2hR2trDIFmax2−2hT2trDIFmax2−2rDIFmax2yT2t+hR4t−2hR2thT2t−2hR2tyT2t+hT4t+2hT2tyT2t+yT4t.In the SAS, since we havehT=hR, (18) and (19) can be simplified as (20) and (21), respectively:(20)μ2t=4rDIFmax2yTt−4yT3t,(21)μ3t=rDIFmax4−4hR2trDIFmax2−2rDIFmax2yT2t+2hR4t−2hR2thT2+yT4t.Thus, the PDF of the scatterers can be denoted as (22) for UD of the scatters in the abovementioned scenarios:(22)fXDYDxD,yD,t=1πrDmax2t,VPS,rDIFmax4μ1tyDmaxt+μ2t⋅μ1tyDmax2t+μ2tyDmaxt+μ3t+logμ1tyDmax2t+μ2tyDmaxt+μ3t+μ2t+2μ1tyDmaxt/2μ1t8μ13/2t⋅4μ1tμ3t−μ22t−4μ1tyDmint+μ2t⋅μ1tyDmin2t+μ2tyDmint+μ3t−logμ1tyDmin2t+μ2tyDmint+μ3t+μ2t+2μ1tyDmint/2μ1t8μ13/2t⋅4μ1tμ3t−μ22t−1,GS&SAS.


## 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 by(23)ht,τ=hLOSt,τ+hSPEt,τ+hDIFt,τ,where hLOSt,τ is the impulse response of the LOS path, hSPEt,τ is the impulse response of the specular reflection path, and hDIFt,τ 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(24)hLOSt,τ=αLOStgLOStδτ−τLOS,where αLOSt=λGTLOStGRLOSt/4πdLOSt, with the assumption that the propagation loss attenuation factor equals to 2. λ is the wave length, and GTLOSt and GRLOSt are the antenna gains of the transmitter and the receiver at the moment t. τLOSt=dLOSt/c, where dLOSt is the length of the LOS path as in (1) at the moment t. gLOSt in (24) is calculated as(25)gLOSt=ejk0dLOSt+j2π/λ∫0tnRTTt⋅vRTtdt,where k0=2π/λ, vRTt=vRt−vTt, vTt is the velocity vector of the transmitter, vTt=vT_xt,vT_yt,vT_ztT, vRt is the velocity vector of the receiver, and vRt=vR_xt,vR_yt,vR_ztT. The unit directional vector of the LOS path from the transmitter towards the receiver is nRTt=nRT_xt,nRT_yt,nRT_ztT. The impulse response of the specular reflection path is given by(26)hSPEt,τ=Γt⋅αSPEt⋅gSPEt⋅δτ−τSPE,where τSPE=dSPEt/c. αSPEt and gSPEt can be calculated by (27) and (28):(27)αSPEt=λGTSPEtGRSPEt4πdSPEt,(28)gSPEt=ejk0dSPEt+j2π/λ∫0tnTSTt⋅vTt+nRSTt⋅vRtdt,where nTS and nRS 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 isL; then, the impulse response of the scattering path is as(29)hDIFt,τ=∑l=1LαDIF,ltgDIF,ltδτ−τDIF,l.When the radar cross section (RCS)σ is constant for each scatterer, the amplitude attenuation of the lth subpath αlt is calculated by (30) in accord with [21, 35]. In (30), dTD,lt and dDR,lt are the distance between the transmitter and the scatterer, and the distance between the scatterer and the receiver, respectively:(30)αDIF,lt=λGTDIFtGRDIFt4πdTD,ltdDR,ltσ.Thus,gDIF,lt in (29), containing the phase and frequency offset in the impulse response of the diffuse scattering components, can be denoted as (31):(31)gDIF,lt=ejφl+jk0dDIF,lt⋅ej2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt,where nTD,lt=nTD,l_xt,nTD,l_yt,nTD,l_ztT and nRD,lt=nRD,l_xt,nRD,l_yt,nRD,l_ztT are the unit directional vectors towards the lth diffuse scatterer D observed from T and R at the moment t, respectively. Thus, the TF and ACF of the scattering paths are calculated as (32) and (33):(32)HDIFt,f=∑l=1LαDIF,lte−j2πfτDIF,ltejφlt⋅ejk0dDIF,lt⋅ej2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt,(33)RDIFt,Δt;f,Δf=EHDIF∗t,f⋅HDIFt+Δt,f+Δf=E∑l=1LλGT,lDIFtGR,lDIFt4πdTD,ldRD,lej2πfτDIF,lte−jφlte−jk0dDIF,lte−j2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt∑l=1LλGT,lDIFt+ΔtGR,lDIFt+Δt4πdTD,ldDR,le−j2πf+ΔfτDIF,lejφltejφlt+Δtejk0dDIF,lt+Δt⋅ej2π/λ∫0t+ΔtnTD,lTt⋅vTt+nRD,lTt⋅vRtdt.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):(34)HDIFt=∑l=1LαDIF,lte−j2πfcτDIF,ltejφlt⋅ejk0dDIF,lt⋅ej2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt,(35)RDIFt,Δt=EHDIF∗t⋅HDIFt+Δt=E∑l=1LλGT,lDIFtGR,lDIFt4πdTD,ldRD,le​jφlte−jk0dDIF,lte−j2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt⋅∑l=1LλGT,lDIFt+ΔtGR,lDIFt+Δt4πdTD,ldDR,le​jφlt+Δtejk0dDIF,lt+Δt⋅ej2π/λ∫0t+ΔtnTD,lTt⋅vTt+nRD,lTt⋅vRtdt.The time-variant PSD of the scattering path can be obtained through Fourier transform with (33) as (36), where u denotes Δt:(36)PDSt,f=∫RDIFt,u⋅exp−j2πfudu.


### 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 Figure2, at the moment t, the length of the scattering path rDIFt can be denoted as (37), where rDt is the radius of the scattering area in VPS at the moment t:(37)rDIFt=hR2t+rD2t+hT2t+rD2t.Thus, we have(38)rDt=ξt2rDIFt.At this moment, the Jacobian [36] can be calculated as(39)J=drDdrDIF=ξ2rDIF2−14rDIFξ⋅hR+hT−rDIFhR−hT+rDIFhR+hT+rDIF+hR+hT−rDIFhT−hR+rDIFhR+hT+rDIF−hR−hT+rDIFhT−hR+rDIFhR+hT+rDIF+hR+hT−rDIFhR−hT+rDIFhT−hR+rDIF,where(40)ξ=−hR+hT−rDIFhR−hT+rDIF⋅hT−hR+rDIFhR+hT+rDIF.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:(41)fRDΦDrD,ϕD=rDπrDmax2,rD∈0,rDmax,ϕD∈0,2π.Thus, the CDF of path delayτDIF can be calculated with (42):(42)PTDIF≤τDIF=PSRDIF≤SrDIF=∫02πdϕ∫0rDrdrπrDmax2=rD2rDmax2=1/2rDIF−hR+hT−rDIFhR−hT+rDIFhT−hR+rDIFhR+hT+rDIF21/2rDIFmax−hR+hT−rDIFmaxhR−hT+rDIFmaxhT−hR+rDIFmaxhR+hT+rDIFmax2=rDIFmax2rDIF2⋅hR+hT−rDIFhR−hT+rDIFhT−hR+rDIFhR+hT+rDIFhR+hT−rDIFmaxhR−hT+rDIFmaxhT−hR+cτmaxhR+hT+rDIFmax=rDIFmax2cτDIF2⋅hR+hT−cτDIFhR−hT+cτDIFhT−hR+cτDIFhR+hT+cτDIFhR+hT−rDIFmaxhR−hT+rDIFmaxhT−hR+rDIFmaxhR+hT+rDIFmax.For the most general scenario, i.e., GS, the CDF ofτDIF is calculated with (43)–(48) and can be simplified with hT=hR for SAS:(43)FTDIFτDIF=PTDIF≤τDIF=∬TDIF≤τDIFfXDYDxD,yDdxDdyD=1SE∬TDIF≤τDIFdxDdyD=STτDIFSE=1SEcτDIF⋅4μ1yD_upper+μ2μ1yD_upper2+μ2yD_upper+μ3+logμ1yD_upper2+μ2yD_upper+μ3+μ2+2μ1yD_upper/2μ14μ1μ3−μ228μ13/2−4μ1yD_lower+μ2μ1yD_lower2+μ2yD_lower+μ3−logμ1yD_lower2+μ2yD_lower+μ3+μ2+2μ1yD_lower/2μ14μ1μ3−μ228μ13/2,(44)yD_upper=hT2yT−hR2yT+yT3−c2τDIF2yT+cτDIF−c2τDIF2+hR2−2hRhT+hT2+yT2⋅−c2τDIF2+hR2+2hRhT+hT2+yT22yT2−2c2τDIF2,(45)yD_lower=hT2yT−hR2yT+yT3−c2τDIF2yT−cτDIF−c2τDIF2+hR2−2hRhT+hT2+yT2⋅−c2τDIF2+hR2+2hRhT+hT2+yT22yT2−2c2τDIF2,(46)μ1=−4c2τDIF2+4yT2,(47)μ2=4c2τDIF2yT+4hR2yT−4hT2yT−4yT3,(48)μ3=c4τDIF4−2c2hR2τDIF2−2c2hT2τDIF2−2c2τDIF2yT2+hR4−2hR2hT2−2hR2yT2+hT4+2hT2yT2+yT4.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. The distributions of the AOAs have been presented in [37].


## 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(23)ht,τ=hLOSt,τ+hSPEt,τ+hDIFt,τ,where hLOSt,τ is the impulse response of the LOS path, hSPEt,τ is the impulse response of the specular reflection path, and hDIFt,τ 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(24)hLOSt,τ=αLOStgLOStδτ−τLOS,where αLOSt=λGTLOStGRLOSt/4πdLOSt, with the assumption that the propagation loss attenuation factor equals to 2. λ is the wave length, and GTLOSt and GRLOSt are the antenna gains of the transmitter and the receiver at the moment t. τLOSt=dLOSt/c, where dLOSt is the length of the LOS path as in (1) at the moment t. gLOSt in (24) is calculated as(25)gLOSt=ejk0dLOSt+j2π/λ∫0tnRTTt⋅vRTtdt,where k0=2π/λ, vRTt=vRt−vTt, vTt is the velocity vector of the transmitter, vTt=vT_xt,vT_yt,vT_ztT, vRt is the velocity vector of the receiver, and vRt=vR_xt,vR_yt,vR_ztT. The unit directional vector of the LOS path from the transmitter towards the receiver is nRTt=nRT_xt,nRT_yt,nRT_ztT. The impulse response of the specular reflection path is given by(26)hSPEt,τ=Γt⋅αSPEt⋅gSPEt⋅δτ−τSPE,where τSPE=dSPEt/c. αSPEt and gSPEt can be calculated by (27) and (28):(27)αSPEt=λGTSPEtGRSPEt4πdSPEt,(28)gSPEt=ejk0dSPEt+j2π/λ∫0tnTSTt⋅vTt+nRSTt⋅vRtdt,where nTS and nRS 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 isL; then, the impulse response of the scattering path is as(29)hDIFt,τ=∑l=1LαDIF,ltgDIF,ltδτ−τDIF,l.When the radar cross section (RCS)σ is constant for each scatterer, the amplitude attenuation of the lth subpath αlt is calculated by (30) in accord with [21, 35]. In (30), dTD,lt and dDR,lt are the distance between the transmitter and the scatterer, and the distance between the scatterer and the receiver, respectively:(30)αDIF,lt=λGTDIFtGRDIFt4πdTD,ltdDR,ltσ.Thus,gDIF,lt in (29), containing the phase and frequency offset in the impulse response of the diffuse scattering components, can be denoted as (31):(31)gDIF,lt=ejφl+jk0dDIF,lt⋅ej2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt,where nTD,lt=nTD,l_xt,nTD,l_yt,nTD,l_ztT and nRD,lt=nRD,l_xt,nRD,l_yt,nRD,l_ztT are the unit directional vectors towards the lth diffuse scatterer D observed from T and R at the moment t, respectively. Thus, the TF and ACF of the scattering paths are calculated as (32) and (33):(32)HDIFt,f=∑l=1LαDIF,lte−j2πfτDIF,ltejφlt⋅ejk0dDIF,lt⋅ej2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt,(33)RDIFt,Δt;f,Δf=EHDIF∗t,f⋅HDIFt+Δt,f+Δf=E∑l=1LλGT,lDIFtGR,lDIFt4πdTD,ldRD,lej2πfτDIF,lte−jφlte−jk0dDIF,lte−j2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt∑l=1LλGT,lDIFt+ΔtGR,lDIFt+Δt4πdTD,ldDR,le−j2πf+ΔfτDIF,lejφltejφlt+Δtejk0dDIF,lt+Δt⋅ej2π/λ∫0t+ΔtnTD,lTt⋅vTt+nRD,lTt⋅vRtdt.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):(34)HDIFt=∑l=1LαDIF,lte−j2πfcτDIF,ltejφlt⋅ejk0dDIF,lt⋅ej2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt,(35)RDIFt,Δt=EHDIF∗t⋅HDIFt+Δt=E∑l=1LλGT,lDIFtGR,lDIFt4πdTD,ldRD,le​jφlte−jk0dDIF,lte−j2π/λ∫0tnTD,lTt⋅vTt+nRD,lTt⋅vRtdt⋅∑l=1LλGT,lDIFt+ΔtGR,lDIFt+Δt4πdTD,ldDR,le​jφlt+Δtejk0dDIF,lt+Δt⋅ej2π/λ∫0t+ΔtnTD,lTt⋅vTt+nRD,lTt⋅vRtdt.The time-variant PSD of the scattering path can be obtained through Fourier transform with (33) as (36), where u denotes Δt:(36)PDSt,f=∫RDIFt,u⋅exp−j2πfudu.


## 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 Figure2, at the moment t, the length of the scattering path rDIFt can be denoted as (37), where rDt is the radius of the scattering area in VPS at the moment t:(37)rDIFt=hR2t+rD2t+hT2t+rD2t.Thus, we have(38)rDt=ξt2rDIFt.At this moment, the Jacobian [36] can be calculated as(39)J=drDdrDIF=ξ2rDIF2−14rDIFξ⋅hR+hT−rDIFhR−hT+rDIFhR+hT+rDIF+hR+hT−rDIFhT−hR+rDIFhR+hT+rDIF−hR−hT+rDIFhT−hR+rDIFhR+hT+rDIF+hR+hT−rDIFhR−hT+rDIFhT−hR+rDIF,where(40)ξ=−hR+hT−rDIFhR−hT+rDIF⋅hT−hR+rDIFhR+hT+rDIF.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:(41)fRDΦDrD,ϕD=rDπrDmax2,rD∈0,rDmax,ϕD∈0,2π.Thus, the CDF of path delayτDIF can be calculated with (42):(42)PTDIF≤τDIF=PSRDIF≤SrDIF=∫02πdϕ∫0rDrdrπrDmax2=rD2rDmax2=1/2rDIF−hR+hT−rDIFhR−hT+rDIFhT−hR+rDIFhR+hT+rDIF21/2rDIFmax−hR+hT−rDIFmaxhR−hT+rDIFmaxhT−hR+rDIFmaxhR+hT+rDIFmax2=rDIFmax2rDIF2⋅hR+hT−rDIFhR−hT+rDIFhT−hR+rDIFhR+hT+rDIFhR+hT−rDIFmaxhR−hT+rDIFmaxhT−hR+cτmaxhR+hT+rDIFmax=rDIFmax2cτDIF2⋅hR+hT−cτDIFhR−hT+cτDIFhT−hR+cτDIFhR+hT+cτDIFhR+hT−rDIFmaxhR−hT+rDIFmaxhT−hR+rDIFmaxhR+hT+rDIFmax.For the most general scenario, i.e., GS, the CDF ofτDIF is calculated with (43)–(48) and can be simplified with hT=hR for SAS:(43)FTDIFτDIF=PTDIF≤τDIF=∬TDIF≤τDIFfXDYDxD,yDdxDdyD=1SE∬TDIF≤τDIFdxDdyD=STτDIFSE=1SEcτDIF⋅4μ1yD_upper+μ2μ1yD_upper2+μ2yD_upper+μ3+logμ1yD_upper2+μ2yD_upper+μ3+μ2+2μ1yD_upper/2μ14μ1μ3−μ228μ13/2−4μ1yD_lower+μ2μ1yD_lower2+μ2yD_lower+μ3−logμ1yD_lower2+μ2yD_lower+μ3+μ2+2μ1yD_lower/2μ14μ1μ3−μ228μ13/2,(44)yD_upper=hT2yT−hR2yT+yT3−c2τDIF2yT+cτDIF−c2τDIF2+hR2−2hRhT+hT2+yT2⋅−c2τDIF2+hR2+2hRhT+hT2+yT22yT2−2c2τDIF2,(45)yD_lower=hT2yT−hR2yT+yT3−c2τDIF2yT−cτDIF−c2τDIF2+hR2−2hRhT+hT2+yT2⋅−c2τDIF2+hR2+2hRhT+hT2+yT22yT2−2c2τDIF2,(46)μ1=−4c2τDIF2+4yT2,(47)μ2=4c2τDIF2yT+4hR2yT−4hT2yT−4yT3,(48)μ3=c4τDIF4−2c2hR2τDIF2−2c2hT2τDIF2−2c2τDIF2yT2+hR4−2hR2hT2−2hR2yT2+hT4+2hT2yT2+yT4.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. 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 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. 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 t=0 s, we assume that the transmitter and the receiver fly towards each other with the velocity of the transmitter vT=0,68,0T m/s in Cartesian coordinates and the velocity of the receiver vR=0,−68,0T m/s. Thus, at the moment t=5 s, the scenario becomes VPS described in Section 3. At the moment t=5 s, we assume that vT=0,−68,0T m/s and vR=0,68,−61T m/s. Finally, at the moment t=10 s, the aircrafts become a SAS and still with the same velocities. The dielectric coefficients are set in accord with [33, 34].Table 2
Simulation parameters.
ParameterValueInitial altitude of the transmitter305 mInitial altitude of the receiver610 mInitial LOS distance between the transmitter and the receiver745 mτmax3.57τLOSCentral frequency1GHzVelocity vector of the transmittert=0 s 0,68,0T m/s, t=5 s, and t=10 s 0,−68,0T m/sVelocity vector of the receivert=0 s 0,−68,0T m/s, t=5 s, and t=10 s 0,68,−61T m/s


### 4.2. Simulation Results of CIR, TF, ACF, and PSD

In our simulations of the CIRs, the tapped delay lines model is adopted. Figure4 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.Figure 4
The CIRs att=0 s and t=5 s (case 1).It can be seen from Figure4 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 duringt∈0,0.12s is as shown in Figure 5, which depicts the time-variant feature of HDIFt,f.Figure 5
The time-frequency channel transfer function of the scattering path.The normalized ACFs of the scattering path at the moments oft=0s and t=5s are shown in Figure 6.Figure 6
The normalized ACFs of the scattering path.From Figure6, we can observe that the main lobe of the normalized ACF at t=0s is much wider than that at t=5s. 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 HDIFt,f.The normalized Doppler PSDs at the momentst=0 s and t=5 s are shown in Figure 7.Figure 7
The normalized Doppler PSD.Figure7 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.


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

The 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.Figure 8
The CDFs ofτDIF.From Figure8, we can see that CDF of τDIF 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 t=0 s and t=10 s.The 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.Figure 9
The CDFs ofαDIF with ϵ=3.57.From Figure9, 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 a1=1.858E+06,b1=3.484E−07,c1=2.308E−08,a2=1.361E+06,b2=3.861E−07,c2=3.81E−08,a3=1.011E+06,b3=4.423E−07,c3=7.004E−08,a4=6.59E+05,b4=5.384E−07,c4=1.221E−07, a5=1.174E+06,b5=8.48E−07,c5=2.817E−07,a6=−8.275E+05,b6=8.915E−07,c6=2.328E−07, and a7=1.581E+05,b7=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].Figure 10
The PDF ofαDIF and its fitted function (with t = 0 s, GS, and ϵ=3.57).(49)fADIFαDIF=∑k=17ake−αDIF−bk2/ck2.The PDF ofαDIF 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 αDIF.Figure 11
The PDF ofαDIF and its fitted function (with t = 5s, VPS, and ϵ=3.57).(50)fADIFαDIF=1.253∗107e−1.089∗106αDIF+3881e1.092∗106αDIF.The simulated PDF ofαDIF in SAS, i.e., at the moment t=10 s, is presented in Figure 12 and is fitted with the 4th order Gaussian function, as (51), where a1=3.033E+06,b1=2.179E−07,c1=1.003E−08,a2=1.658E+06,b2=2.46E−07,c2=2.846E−08, a3=8.147E+05,b3=2.934E−07,c3=6.558E−08,anda4=1.99E+06,b4=1.17E−07,c4=5.041E−07. Compared with the GS, a simpler expression of the PDF of αDIF is obtained with best goodness of fit.Figure 12
The PDF ofαDIF and its fitted function (with t = 10 s, SAS, and ϵ=3.57).(51)fADIFαDIF=∑k=14ake−αDIF−bk2/ck2.


## 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 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. 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 t=0 s, we assume that the transmitter and the receiver fly towards each other with the velocity of the transmitter vT=0,68,0T m/s in Cartesian coordinates and the velocity of the receiver vR=0,−68,0T m/s. Thus, at the moment t=5 s, the scenario becomes VPS described in Section 3. At the moment t=5 s, we assume that vT=0,−68,0T m/s and vR=0,68,−61T m/s. Finally, at the moment t=10 s, the aircrafts become a SAS and still with the same velocities. The dielectric coefficients are set in accord with [33, 34].Table 2
Simulation parameters.
ParameterValueInitial altitude of the transmitter305 mInitial altitude of the receiver610 mInitial LOS distance between the transmitter and the receiver745 mτmax3.57τLOSCentral frequency1GHzVelocity vector of the transmittert=0 s 0,68,0T m/s, t=5 s, and t=10 s 0,−68,0T m/sVelocity vector of the receivert=0 s 0,−68,0T m/s, t=5 s, and t=10 s 0,68,−61T m/s


## 4.2. Simulation Results of CIR, TF, ACF, and PSD

In our simulations of the CIRs, the tapped delay lines model is adopted. Figure4 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.Figure 4
The CIRs att=0 s and t=5 s (case 1).It can be seen from Figure4 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 duringt∈0,0.12s is as shown in Figure 5, which depicts the time-variant feature of HDIFt,f.Figure 5
The time-frequency channel transfer function of the scattering path.The normalized ACFs of the scattering path at the moments oft=0s and t=5s are shown in Figure 6.Figure 6
The normalized ACFs of the scattering path.From Figure6, we can observe that the main lobe of the normalized ACF at t=0s is much wider than that at t=5s. 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 HDIFt,f.The normalized Doppler PSDs at the momentst=0 s and t=5 s are shown in Figure 7.Figure 7
The normalized Doppler PSD.Figure7 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.


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

The 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.Figure 8
The CDFs ofτDIF.From Figure8, we can see that CDF of τDIF 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 t=0 s and t=10 s.The 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.Figure 9
The CDFs ofαDIF with ϵ=3.57.From Figure9, 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 a1=1.858E+06,b1=3.484E−07,c1=2.308E−08,a2=1.361E+06,b2=3.861E−07,c2=3.81E−08,a3=1.011E+06,b3=4.423E−07,c3=7.004E−08,a4=6.59E+05,b4=5.384E−07,c4=1.221E−07, a5=1.174E+06,b5=8.48E−07,c5=2.817E−07,a6=−8.275E+05,b6=8.915E−07,c6=2.328E−07, and a7=1.581E+05,b7=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].Figure 10
The PDF ofαDIF and its fitted function (with t = 0 s, GS, and ϵ=3.57).(49)fADIFαDIF=∑k=17ake−αDIF−bk2/ck2.The PDF ofαDIF 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 αDIF.Figure 11
The PDF ofαDIF and its fitted function (with t = 5s, VPS, and ϵ=3.57).(50)fADIFαDIF=1.253∗107e−1.089∗106αDIF+3881e1.092∗106αDIF.The simulated PDF ofαDIF in SAS, i.e., at the moment t=10 s, is presented in Figure 12 and is fitted with the 4th order Gaussian function, as (51), where a1=3.033E+06,b1=2.179E−07,c1=1.003E−08,a2=1.658E+06,b2=2.46E−07,c2=2.846E−08, a3=8.147E+05,b3=2.934E−07,c3=6.558E−08,anda4=1.99E+06,b4=1.17E−07,c4=5.041E−07. Compared with the GS, a simpler expression of the PDF of αDIF is obtained with best goodness of fit.Figure 12
The PDF ofαDIF and its fitted function (with t = 10 s, SAS, and ϵ=3.57).(51)fADIFαDIF=∑k=14ake−αDIF−bk2/ck2.


## 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].



---

*Source: 9267907-2022-01-13.xml*