Analytical Non-Stationary Satellite to Aircraft Channel Modeling over Open Area Based on Regular Shaped Geometry-Based Stochastic Model

: Channel modeling is crucial to the development and evaluation of modern wireless communication systems including satellite communication system, since there might be critical safty-of-life applications. Also, the channel model is of great importance to the performance evaluation of mobile communication systems. In recent years, encouraged by the widely application of unmanned aerial vehicles, the research on channel modeling for aerial and aeronautical communications attract lots of interests. In the published articles, stationary and non-stationary channel models have been developed for air-to-ground communications based on regular shaped geometry-based stochastic model (RS-GBSM). The modeling of air-to-air or satellite-to-aircraft (S2A) communication is still quite simple or completely lacking. For obtaining more precise model of S2A channel, this paper presents an analytical non-stationary S2A channel mode based on RS-GBSM with considerations on line-of-sight path, specular reﬂection path, and ground scattering path. Analytical expressions of the channel impulse responses, the transfer functions, the auto-correlation functions, and the Doppler power spectrum density based on 3-path model are derived and simulated. Also, the distributions of the path antennation, the path delay, and the normalized Doppler shift based on uniform distribution of the scatterers are derived, simulated and ﬁtted.


Introduction
The accuracy of performance evaluation on satellite-to-aircraft (S2A) communication is important since there might be critical safety-of-life applications. Moreover, the accuracy of system performance evaluation is to some extend determined by precise modeling of the channel. Existing techniques for modeling of the S2A channel are not sufficiently precise [1], although some literatures such as [1][2][3][4] have proposed methodologies on S2A channel modeling or have presented measurements to improve the channel models. Fortunately, there are some preliminary studies on the channel modeling of air-to-ground (A2G) communications, air-to-air (A2A) communications, and vehicle-to-vehicle (V2V) communications that are worth learning. For example, the methods of developing regular shaped geometry-based stochastic models (RS-GBSMs) were presented in [5][6][7], and non-stationary characteristics of mobile channels were discussed in [8]. With the methods of developing non-stationary RS-GBSM and introducing the ground scattering that must be considered in S2A channel modeling, we are able to develop and present an analytical non-stationary RS-GBSM based S2A model in this paper. exponential function. δ (·) is the unit impulse function. P(X) is the probability of event X happening. F(·) is the CDF. f (·) is the PDF. h (·) is the CIR. H (·) is the TF. R (·) is the ACF.

The Geometrical Structure
The presented structure is as shown in Figure 1, with Figure 1a presenting the system structure with spherical earth surface, and Figure 1b illustrating the paths and scattering area. Major parameters used in Figure 1 and further formulations are also listed in Table 1.   Time variant distance between the transmitter and the receiver (length of the LOS path). d SPE (t) The time variant length of the specular reflection path d TO (t) The time variant distance between the transmitter T and the specular reflection point O d OR (t) The time variant distance between the specular reflection point O and the receiver R d DIF,l (t) The time variant length of the lth scattering sub-path.
The antenna gain of the receiver at the moment t along the LOS path G LOS R (t) The antenna gain of the receiver at the moment t along the LOS path G SPE T (t) The antenna gain of the receiver at the moment t along the specular reflection path G SPE R (t) The antenna gain of the receiver at the moment t along the specular reflection path G DIF T,l (t) The antenna gain of the receiver at the moment t along the specular reflection path G DIF R,l (t) The antenna gain of the receiver at the moment t along the specular reflection path h R ' The altitude of the receiver h R (t) The time variant altitude of the receiver h T ' The altitude of the transimitter h T (t) The time variant altitude of the transimitter h T The height of the transimitter h R The height (the Z coordinate of R) h R (t) The height (the Z coordinate of R) h T The height (the Z coordinate of T) h T (t) The time variant height (the Z coordinate of T) h (t, τ) CIR of the channel h LOS (t, τ) CIR of the LOS path h SPE (t, τ) CIR of the specular reflection path h DIF (t, τ) CIR of the scattering path l l (t) The length of the half major aix of the scattering ellipse l s (t) The length of the half minor aix of the scattering ellipse The unity directional vector of the LOS path from the transmitter towards the receiver n TO The unit directional vector towards O observed from T n RO The unit directional vector towards O observed from R n TD,l The unit directional vector towards the lth scatterer observed from the transmitter T n RD,l The unit directional vector towards the lth scatterer observed from the receiver R O' The center of the earth O The origin of the XOY plane, also the specular reflection point r e The radius of the earth R The receiver (aircraft) The time variant area of the scattering ellipse t Time T The The time variant 3D coordinate of the transmitter The time variant 3D coordinate of the receiver The time variant 3D coordinate of the scatterer (x D,l (t) , y D,l (t) , 0 (t)) The time variant 3D coordinate of the lth scatterer α 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 ∆ f Frequency interval ε e r Relative dielectric permeittivity γ 1 The angle between O'R and Z axis γ 2 The angle between O'T and Z axis Time variant reflection for vertical polarization Γ // (t) Time variant reflection coefficient for horizontal polarization γ DIF Normalized Doppler shift of the scattering path Φ The geocentric angle between T (the satellite) and R (the aircraft) The geocentric angle between T (the satellite) and R (the aircraft) φ l Additional phase shift of the lth scattering path σ Radar cross section (RCS) τ Delay τ LOS (t) Time variant path delay of the LOS path at the moment t τ SPE (t) Time variant path delay of the specular reflection path at the moment t τ DIF (t) Time variant path delay of the scattering path at the moment t τ DIF,l (t) Time variant path delay of the lth scattering path at the moment t θ Reflection angle θ (t) Time variant reflection angle In Figure 1, T is the satellite transmitter, R is the aircraft receiver. The altitudes of R and T are h R ' and h T ' respectively. O' is the center of the earth, and O is the ground specular reflection point. In this paper, O is determined with C.Wagner algorithm [22], and is set as the origin of the Cartesian coordinates. θ = 1 2 ROT is the reflection angle. The YOZ plane is determined with R, O, and T, and the Z axis is set alone O O. Y axis is set to be perpendicular to Z axis and directing to T. h T and h R are the heights and also the Z axis coordinates of T and R. Φ = γ 1 + γ 2 = RO'T. According to the concerned literatures, such as [2,23], the earth surface should be modeled as the spherical surface. Since the atmospheric refraction should be considered, effective earth radius r e is adopted, as shown in Figure 1. In this paper, the longitudes, latitudes, and altitudes of the satellite and the aircraft are assumed to be known. And the satellite is assumed to be a Geosynchronous Earth Orbit (GSO) satellite. Similar to [2], the effects of tropospheric are not considered in this paper since ground scattering is what we focus on and the effects of tropospheric is quite stable in the scenario we discuss according to [24].
In Figure 1, since the transmitter and the receiver might moving during communication, the distances and the AOAs are time variant. Thus, in most of the situations the channel is non-stationary. According to the common method of developing aeronautical channel models, such as presented in [2], three paths, including the LOS path, the reflection path, and the diffuse scattering path, are considered in this paper. Based on Figure 1, the lengths of the LOS path and the specular reflection path can be calculated by Equations (1) and (2).
As shown in Figure 1, in this paper, for simplifying the modeling, we assume that the scattering happens on a horizontal surface intersects with the earth surface at the specular reflection point O, thus the coordinates of the transmitter, the receiver, the scatterer are known as Since the scattering path contains sub-paths, when we discuss the sub-paths, the coordinate of the lth scatterer is denoted as (x D,l (t) , y D,l (t) , 0), thus the length of the scattering sub-path can be calculated as Equation (3).
In this paper, the scatterers are assumed to be distributed in an elliptical region determined by the maximum delay ellipsoid with T and R as its foci, and truncated by the XOY plane with the specular reflection point as the center, as shown in Figure 1. At the moment t, the range of the scattering area is constrained by Equation (4),where τ max is the maximum path delay compared to the LOS path, and c is the speed of light.
For simplification, we take this scattering region constrained by Equation (4) as an ellipse with half major axis length l l and half minor axis length l s , where l l and l s are calculated through resolving Equation (4) with conditions x D (t) = 0 and y D (t) = 0, respectively, and as shown in Equations (6) and (7), where Thus the area of the scattering region is calculated as in Equation (8).
It's to be mentioned that, the error of calculating the region area through this approximation, i.e., taking the scattering region as an ellipse, is under 0.1% in the presented scenarios. Calculation examples can be seen in Appendix A.
We assume that the scatterers are obey regional UD, thus we have the joint PDF of the time variant X coordinate x D (t), and the time variant Y coordinate y D (t) of the scatterer at the moment t, as in Equation (9).

Channel Impulse Response
As mentioned above, three paths model is adopted in this paper, thus the CIR of the presented channel model can be given by Equation (10).
where h LOS (t, τ) is the impulse response of the LOS path with t as the observed moment, and τ as the corresponding time delay, 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. Considering the large scale propagation attenuation, propagation delay and Doppler shift, the impulse response of the LOS path is denoted is the path attenuation of the LOS path, 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 corresponding g LOS (t) in Equation (11), is the term for phase shifts introduced by LOS transmitting and Doppler shift, and is calculated as in Equation (12), T is the velocity of the receiver. The unity 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 .
The impulse response of the specular reflection path is given by Equation (13) , where τ SPE = d SPE (t)/c. α SPE (t) and g SPE (t) can be calculated by Equations (14) and (15).
In Equations (14) and (15), n TO and n RO are the unit directional vector towards O observed from T and R, respectively. And Γ (t) can be calculated as Fresnel Equations (16) and (17), with vertical and horizontal polarization, respectively,where ε e r (t) is the relative dielectric permittivity of the reflection surface, which is a random process since the reflector and the scatterers may be changed with the movements of the transmitter and the receiver, the reflecting and scattering areas are for example, the sea, earth surface, etc. [25][26][27].
In the simulation of this paper, ε e r (t) is calculated as in [27] with the assumption that the scattering surface is sea surface.
Suppose that the number of the scatterers is L, then the impulse response of the scattering path is as Equation (18).
When the radar cross section (RCS) σ is constant for each scatterer, the amplitude attenuation of the l th subpath α l (t) is calculated by Equation (19) in accord with [2,28]. In Equation (19), 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.
According to the radar equation and the theory presented in [2], when RCS is not constant in the diffuse scattering area, we can denote it as σ (x, y). Thus α DIF,l (t) is calculated as in (20) with the assumption that the shape of each scatterer is a circle C l with radius R l .
where σ x and σ y are the root mean square of the slop of x and y, respectively, b xy is the correlation coefficient of x and y, and k is calculated as in Equation (23).
Thus the term g DIF,l (t) in Equation (18), containing the phase and frequency shifts in the impulse response of the scattering components, can be denoted as Equation (24), where n TD,l = n TD,l_x , n TD,l_y , n TD,l_z T and n RD,l = n RD,l_x , n RD,l_y , n RD,l_z T are the unit directional vectors towards the lth diffuse scatterer D observed from T and R, separately. ϕ l ∼ U [−π, π] is an additional random phase for the lth scattering subpath.

Time Variant Transfer Function
The time variant transfer function H (t, f ) can be obtained by Fourier transform of impulse response h (t, τ) with respect to the variable τ , as shown in Equation (25). With H (t, f ), time variant random spectrum of the channel can be observed.
The terms H LOS (t, f ), H SPE (t, f ), and H DIF (t, f ) denote the time-variant transfer functions of the LOS component, specular reflection component, and scattering components, respectively. Thus Equations (26)-(28) can be obtained.
Equations (26)- (28) can be used to derive time-frequency auto-correlation functions as in Section 3.3.

Time-Frequency ACFs
With the assumption that the LOS path, the specular reflection path, and the diffuse scattering paths are uncorrelated with each other, the ACFs can be derived through the time variant transfer functions. With time difference denoted as ∆t and frequency difference denoted as ∆ f , the time-frequency correlation function can be expressed as Equation (29) With the assumption that the reflection coefficient Γ, delay τ, distance d, additional phase shift ϕ, antenna gain G, velocity v, and etc. are time independent, the scattering channel satisfies the feature of wide sense stationary uncorrelated scattering (WSSUS) and the auto-correlation functions of each paths can be simplified as Equations (33)- (35) .
With time variant ACFs, different coherence time can be obtained, which is essential parameter for system design. Through Fourier transform [31] as in Equation (36) with Equations (33)- (35), for the WSSUS channel, the Doppler delay spectrum of each path can be obtained, where R (u, ∆ f ) is R LOS (∆t, ∆ f ), R SPE (∆t, ∆ f ), or R DIF (∆t, ∆ f ), u denotes ∆t. If more precise information of the channel is desired, the time variant features should not be ignored. Thus other method towards the frequency spectrum features such as Wigner Ville spectrum should be adopted, for example, as discussed in [13].

Distributions of Path Delay, Path Attenuation, Normalized Doppler Shift of the Scattering Path
In this section, we investigate the time varying distributions of path delay, path attenuation, and normalized Doppler shift of the diffuse scattering path, with the assumption that the scatterers obey UD and the PDF is as Equation (9). The structure of the scattering area is illustrated in Figure 2, where the specular reflection point is the center of the ellipse, l l is the semi-major axis, and l s is the semi-minor axis.

The Distribution of τ DIF (t)
The time variant CDF of the time delay of the scattering path relative to the LOS path can be calculated with Equations (37)-(41) [32], where P (T DIF (t) ≤ τ DIF ) is the probability that the virtual random variable T DIF observed at the moment t is not larger that the random path delay τ DIF observed at the moment t, τ min ≤ τ DIF ≤ τ max , and S T (τ DIF , t) is the area of the elliptical area determined by Equation (38). τ min is the minimum path delay of the scattering path obtained through tangency point between the minimum delay-determined ellipsoid and the earth, τ min = τ SPE .
can be calculated with Equation (39) similar to Equations (6)-(8), as follows, where The details of calculating l l,T (τ DIF , t) and l s,T (τ DIF , t) can be found in Appendix B.
Thus, the PDF f (τ DIF , t) can be obtained with Equation (47), where ∂l l,T (τ DIF ,t) ∂τ DIF is calculated as in (48). where is calculated as in Equation (55), where The calculation of l l,T (τ DIF , t) and l s,T (τ DIF , t) can be refered to Appendix B.

The Distribution of α DIF (t)
The path attenuation of the scattering path can be calculated with Equation (61) according to Radar function [2,28], for any sub-path, it can be denoted as α DIF (t) = q (x D (t) , y D (t)) in Equation (61), where the scatter is assumed to be circle with R l as its radius. The PDF of the scattering path attenuation α DIF (t) can be obtained with Equation (62) according to the theory of random process [32], where S (t) is calculated with Equation (8), F A DIF (α DIF , t) is the CDF of α DIF (t), and A DIF is the virtual random variable at the moment t calculated with Equation (61).

The Distribution of γ DIF (t)
We define the normalized Doppler shift as follows.
It can be observed from Equation (63) that γ DIF is a function of x D (t) and y D (t), thus the PDF of γ DIF (t) can be calculated as Equation (64), where F γ DIF (γ DIF , t) is the time variant CDF of γ DIF (t), and Γ DIF (t) is the virtual random variable calculated with Equation (63).
Since it is very hard to be obtain closed form expressions from Equations (62) and (64), we perform simulation and distribution fitting to get them.

Simulation Assumptions
According to similar scenarios or assumptions used in [1,4,14,[33][34][35], the simulation parameters are set and listed in Table 2. The trajectories of the receiver are shown in Figure 3 for the simulation of the CIRs. The simulation includes three scenarios(states), i.e., rising, cruising, and descending. in the simulation of CIR, TF, ACF, and PSD, we assume that the initial positions of these three scenarios are the same, i.e., the position of R in Figure 3 corresponds to the moment t = 0s in each scenario, then the aircraft moves along the trajectories corresponding to each scenario along the directions of the arrows. The moments we simulate are collected at different moments during the moving of the aircraft. The dielectric coefficients are set in accord with [26,27].
Cruising Figure 3. The Trajectories of the Aircraft Table 2. Simulation Parameters.

Parameter Value
Altitude of the transmitter 36,000 km Initial altitude of the receiver (Rising, t = 0 s) 300 m Initial altitude of the receiver (Cruising, t = 0 s) 11 km Initial altitude of the receiver (Descending, t = 0 s) 11 km Initial horizontal distance between the transmitter and the receiver 400 km τ max τ SPE + 0.7µs − τ LOS Central frequency 2 GHz Ground velocity vector of the receiver (Rising)

Simulations of CIR, TF, ACF, and Normalized Doppler PSD
In our simulations of the CIRs, tapped delay lines model is adopted in this paper with ∆t =0.1 µs. Figure 4 shows a snapshot of the impulse responses of the presented model at the moment t = 1 s of the rising stage. Six scattering taps are illustrated in this figure.
As shown in Figure 4, the LOS path, which is also the path with the minimum delay, can be distinguished from other paths easily, the specular reflection path is adjacent to the scattering paths. Closer scattering components to the specular component would be observed if nearer scatterers to the specular reflection point were distributed.
The time-frequency transfer function of the scattering path during t ∈ (1, 1.01) s of the rising stage is shown in Figure 5 which depicts the time variant feature of the H DIF (t, f ). Although H DIF (t, f ) is only a sample of the spectrum of the random process h DIF (t, τ), we can observe its time variant and random feature.
The normalized ACFs of the scattering path in different flight states and at different moments, i.e., t = 300 s, t = 450 s, t = 650 s, and t = 800 s of the rising state, t = 400 s, t = 500 s, t = 600 s, and t = 700 s of the cruising state, and t = 450 s, t = 600 s, t = 750 s, and t = 900 s of the descending state, are shown and compared in Figures 6-8. From Figure 6, i.e., in the rising state, we can observed that the main lobe of the normalized ACF becomes wider with the increasing of the time, i.e., the distance between the transmitter and the receiver is decreasing. This depicts increasing correlation time periods. We can also observe that the fluctuations of the curves weaken with the increasing of time. These show the non-stationary feature of the channel. In Figure 8 an reverse trend was observed. Quasi-stationary feature can be observed in the cruising state, as shown in Figure 7.     The normalized Doppler PSDs during the time period of t ∈ (0, 1200) s for rising, cruising and descending states are shown in Figure 9, Figure 10, and Figure 11, respectively. These figures show the increases of the maximum normalized Doppler shifts in accord with the movements of the receiver towards the transmitter. As shown in Figure 9, in the rising state, the maximum normalized Doppler becomes larger with the increasing of the time. A reverse trend can be observed in the descending state with Figure 11. Quasi-stationary feature can be observed in the cruising state, as shown in Figure 10.

Simulations of the Distributions of Path Delay, Path Attenuation, and Normalized Doppler Shift of the Scattering Path
In this subsection, simulation results of the distributions of τ DIF are illustrated and compared with the derived CDFs. For the distributions of α DIF and γ DIF , since there are no available closed form expression for the presented scenario, for more convenient application, the fitting functions are presented. For better observation of the symmetries of the Doppler shift, the PDFs are observed instead of CDFs. It's to be noted that, different from the trajectories in Scetion 5.2, the trajectory of the aircraft in descending state is the return route of the trajectory of the rising state, i.e., 1800 s were simulated for both states.

The Distribution of τ DIF
The CDFs of scattering path delay τ DIF in the rising state at the moments of t = 0 s, t = 100 s, and t = 700 s and correspondingly at the moments of t = 1800 s, t = 1700 s, and t = 1100 s, of the descending state are illustrated in Figure 12, as well as the CDFs of τ DIF in cruising state. The simulated CDFs coincide perfectly with the derived CDFs. It's to be mentioned that, for clearer comparison, in Figure 12, τ DIF is illustrated relative to the SPE path, which is also the scattering path with minimum path delay. It can be seen form Figure 12 that, in the rising and the descending states, the CDFs of τ DIF is obviously different at different moments, this depicts the non-stationary feature of the channel. In the rising state, as time increasing, decreasing altitude difference with decreasing distance between the aircraft and the satellite means larger and more circle-like scattering region with the same maximum time delay, i.e., (1.4 µs in this subsection). Thus the CDF of τ DIF tends to be linear. Similar features can be observed with the simulation results at corresponding moments of descending state, since similar geometry structures can be obtained. And the CDFs of τ DIF in the cruising state are quasi-stationary due to almost time invariant geometry of the scatterers, and it is very close to the CDF at the moment of t = 700 s in the rising state.

The Distribution of α DIF
By assuming the distribution of the scatterers is UD, the CDFs of the path attenuation of the scattering path α DIF relative to the minimum scattering path attenuations at the moments t = 0 s, t = 100 s and t = 700 s respectively, of the rising state, which are corresponding to the moments of t = 1800 s, t = 1700 s and t = 1100 s on the descending state, are shown in Figure 13. The CDFs of the path attenuation of the scattering path α DIF relative to the minimum scattering path attenuation at the moments t = 0 s, t = 700 s and t = 1700 s of cruising state are also shown in Figure 13.
From Figure 13, we can observe that, in the rising state, as the distance between the aircraft and the satellite decreasing with t increasing, the path attenuation difference tends to diminish. In the cruising state, quasi-stationary feature can be observed through very close CDFs at different moments. Compared with the features obtained in the rising state, a reverse trend can be observed with the simulation resluts in the descending state. The CDFs of α DIF in rising state can be fitted with two terms Fourier function as shown in (65), with the coefficients listed in Table 3.  The CDFs of α DIF in cruising state can be fitted with one term Fourier function as shown in Equation (66), with the coefficients a α,0 = 0.6214, a α,1 = 0.3768, b α,1 = 0.4866, and w = 2.232.

The Distribution of γ DIF
The PDFs of normalized Doppler shift γ DIF is illustrated in Figures 14-16. In the Rising state, as shown in Figure 14 the PDF of γ DIF becomes narrower with the increasing of observed moments, depicting non-stationary feature of the channel. In Figure 16, a reverse trend can be observed compared with the PDF of γ DIF obtained in the rising state, which also depicts the non-stationary feature of the channel. While as shown in Figure 15 the PDFs of γ DIF in the cruising state are very similar, reflecting quasi-stationary feature of the channel.   The PDFs of γ DIF in the Rising state t = 400 s, t = 600 s and t = 800 s, are fitted with a polynomial as Equation (67). The coefficients are as shown in Table 4. The PDFs of γ in the cruising and descending states can be fitted with the same formula (67), with different coeffcients.

Conclusions
In this paper, an analytical RS-GBSM based S2A channel model for UD of scatterers over an open area is proposed. Based on the presented geometry, the CIRs of the LOS path, the specular reflection path, and the scattering path are formulated. The corresponding time-frequency TFs, ACFs and PSDs are derived and simulated. Simulation scenarios are designed as three states, e.g., the rising state, the cruising state, and the descending state. In each state, several moments are selected to illustrate the time variant features of the channel. Through the simulation results, the non-stationary features of the channel in rising and descending states are evidently observed, while in cruising state, the channel is quasi-stationary. Furthermore, based on an reasonable simplification that the scattering area can be approximated as an ellipse, the closed form expression of path delay distribution is derived, simulation results coincide with the theoretical results perfectly. The distributions of path attenuation, and Doppler shift of the scattering path are also simulated and fitted since no closed form expressions of their CDFs can be derived. Based on the theoretical and simulation results of the distributions of channel coefficients, the non-stationary features in the rising and descending states are observed and analyzed, as well as the quasi-stationary feature in the cruising state. These observations are in accord with the analysis with the channel formulations. Future work will focus on the derivation and analysis based on different distributions of scatterers, for example, von-Mises Fisher (VMF) distribution as adopted in [36]. Since the geometrical structures of channel modeling must base on the realistic environment, further investigation on concerned achievements in geography and physics would be performed. Analytical non-stationary channel models of extended scenarios, e.g., A2A scenario, will be developed as well.
Funding: This research received no external funding.

Conflicts of Interest:
We declare that we have no conflict of interest.

Appendix A. Examples of Area Calculation Error of the Scattering Ellipse
In Table A1 we list some examples of area calculation error of the scattering ellipse, with the presented approximation presented in this paper, the error is under 0.1%.