SATELLITE COMMUNICATION-ORIENTED GEOMETRY-BASED STOCHASTIC CHANNEL MODELING METHOD

Abstract

A satellite communication-oriented geometry-based stochastic channel modeling method includes steps: S1: establishing a satellite channel simulation scenario, and setting scenario layout parameters; S2: initializing the trajectory and speed of a satellite and a receiving end; S3: calculating spatially consistent large-scale parameters, and calculating an effect of rainfall on the large-scale parameters; S4: calculating a path loss, a shadow fading, an atmospheric absorption, and a rainfall attenuation; S5: initializing the central positions of a cluster and a scatterer, and calculating the delay, angle and power of the cluster according to the geometric position information of the transmitting and receiving ends and the scatterer, to generate a channel coefficient; S6: updating the large-scale and small-scale parameters according to the movement of the transmitting and receiving ends and the birth-death process of the cluster, to generate a new channel coefficient; and S7: deriving the statistical characteristics of the channel, and performing simulation analysis.

Description

FIELD OF THE INVENTION

The present invention belongs to the technical field of channel modeling, and in particular, relates to a satellite communication-oriented geometry-based stochastic channel modeling method.

DESCRIPTION OF RELATED ART

In recent years, the fifth generation (5G) mobile communication technology has begun to be commercially deployed around the world, and the research and development of the sixth generation (6G) mobile communication technology is also in full swing. The 6G vision can be summarized as “global-coverage, all-spectra, full-application, full-sense, full-digitalization, and strong security”. It aims to provide higher communication rates, more user connections, and wider network coverage on the basis of 5G, providing users with a more intelligent, secure, and immersive “Everything Connected” experience. In order to achieve deep global coverage, 6G will expand from ground mobile communications to a space-air-ground-sea integrated communication network, of which satellite communication serves an indispensable part. Satellite communication has the advantages of large coverage area, long communication distance, and suitability for a variety of services. It can realize communication with fixed or mobile terminals and is widely used in application scenarios such as navigation, earth observation, and broadcasting. With the rapid development of the space-air-ground-sea integrated network, 6G-oriented satellite communication networks require larger capacity and higher service quality to provide users with ubiquitous connectivity. Typical applications include satellite Internet of Things, satellite-to-cellular service, fixed network backhaul, connected cars, emergency safety communications, aircraft broadcasting, maritime communications, and so on.

Depending on whether the satellite is stationary relative to the earth, satellites can be divided into geostationary earth orbit (GEO) satellites and non-geostationary earth orbit (NGEO) satellites; and according to the altitude of the satellite orbit, satellites can be divided into high earth orbit (HEO) satellites, medium earth orbit (MEO) satellites, and low earth orbit (LEO) satellites. GEO satellites are typical HEO satellites, while LEO and MEO satellites are NGEO satellites. Low-orbit satellites orbit at an altitude of 500-2000 km and have the advantages such as low propagation delay, small size and low signal loss, and are expected to significantly improve the communication speed and energy efficiency. However, in order to achieve global coverage, multiple satellites are required for networking. For example, the StarLink system consists of 1,584 satellites deployed at an altitude of 550 km. MEO satellites have an altitude of 5,000-12,000 km, orbiting at a lower speed, and each has a larger communication range. Satellites at an altitude above 20,000 km are generally HEO satellites. GEO satellites with an altitude of 36,000 km are typical HEO satellites, having an advantage that only three satellites are enough to cover the entire earth, and there is almost no need to track the direction of the satellite. However, the high orbital altitude of the satellite leads to very high free space path loss.

The 6G-oriented satellite communication network shows a trend of evolving towards broadband and Internet. It can be seen from the communication satellites currently in orbit and planned to be launched that future satellite communication systems will adopt millimeter-wave frequency bands such as Ka and Q bands to achieve high-speed communication and broadband digital transmission. At the same time, they will continue to use low-frequency bands that are not easily affected by clouds and rainy weather, to compensate for the shortcomings of the millimeter-wave bands susceptible to weather factors such as rainfall, and provide higher-quality communication services through multiple complementary frequency bands. In terms of orbital altitude, large-scale NGEO broadband Internet constellations are the current focus of global aerospace development, and high-throughput GEO satellites also gain continuous attention thanks to their advantages in coverage, capacity and bandwidth cost. Channel modeling lays a foundation for the design, optimization and adjustment of a communication system, and plays a crucial role in all aspects of a wireless communication system. In view of the development trend that 6G satellite communication networks will adopt multiple frequency bands for communication and satellites at multiple orbital altitudes, a general satellite channel model suitable for different frequency bands and different orbital altitudes needs to be proposed.

SUMMARY OF THE INVENTION

The present invention aims to provide a satellite communication-oriented geometry-based stochastic channel modeling method, which can characterize different channel characteristics of satellite communication systems covering S-band to millimeter-wave band, and is suitable for various satellite trajectories and receiving-end environments to solve the technical problems mentioned in the background technology.

To solve the above-mentioned technical problem, the specific technical solution of the present invention is as follows:

A satellite communication-oriented geometry-based stochastic channel modeling method, wherein the channel modeling method comprises the following steps:

• • Step S1: establishing a satellite communication-oriented geometry-based stochastic channel model and a satellite channel simulation scenario corresponding to the model, and setting scenario layout parameters;
  • Step S2: initializing the trajectory and speed of a satellite and a ground user;
  • Step S3: calculating spatially consistent large-scale parameters, and the effect of rainfall on the large-scale parameters;
  • Step S4: calculating a path loss, a shadow fading, an atmospheric absorption and a rainfall attenuation;
  • Step S5: initializing the central positions of a cluster and a scatterer, and calculating the delay, angle and power of the cluster according to the geometric position information of the transmitting and receiving ends and the scatterer, to generate a channel coefficient;
  • Step S6: updating the large-scale and small-scale parameters according to the movement of the transmitting and receiving ends and the birth-death process of the cluster, to generate a new channel coefficient; and
  • Step S7: deriving the statistical characteristics of the channel, and performing simulation analysis.

Further, the step S1 includes the following specific steps:

• • Step S101: the transmitting end of the model is a satellite, the receiving end is a mobile ground user, the orbital altitude h_sat of the satellite needs to be pre-set according to the requirements of the simulation environment, and the environment of the ground mobile receiving end also needs to be pre-set as a dense city, a city, a suburb, or a rural environment;
  • Step S102: both the transmitting end and the receiving end of the model are equipped with a antenna array, so the spacing and angle parameters of the antenna array need to be initialized, wherein the antenna array is a uniformly distributed linear array; the transmitting end has a total of P antennas, A_p^T denotes the pth antenna, β_A^T denotes an azimuth angle of the transmitting-end antenna array, β_E^T denotes an elevation angle of the transmitting-end antenna array, and a spacing between the antennas is δ^T; and the receiving end has a total of Q antennas, A_q^T denotes the qth antenna, β_A^R denotes an azimuth angle of the receiving-end antenna array, β_E^R denotes an elevation angle of the receiving-end antenna array, and a spacing between the antennas is δ^R.
  • Step S103: the system parameters of the model further comprise a carrier frequency f_c, and a rainfall rate.

Further, the step S2 specifically includes:

• • Step S201: the initial position of the satellite is determined together by the orbital altitude h_sat, elevation θ_sat and azimuth φ_sat; satellites are divided into GEO satellites and NGEO satellites according to whether the satellites are stationary relative to the earth; the GEO satellites are stationary relative to the ground; under the Earth's gravity, the NGEO satellites orbit the earth in an elliptical trajectory; first, an elliptical orbital plane of the satellite is determined by setting the following five parameters: a length a of the semi-major axis of the elliptical orbit; an eccentricity e of the elliptical orbit, wherein if the eccentricity e is set to 0, the elliptical orbit is reduced to a circular orbit; the angle between the orbit and the equatorial plane is an orbital plane inclination ι; the longitude of the point where the orbit crosses the equatorial plane upward is called a right ascension of the ascending node Ω; the direction of the ellipse on the orbital plane is determined by an argument of perigee ω of an angle between the orbital perigee and the ascending node; a true Anomaly υ represents an angle swept by moving along the orbit from the perigee; after determination of the orbital plane, the specific position of the satellite in the elliptical orbit can be determined via υ; the satellite trajectory is determined together by a time-varying right ascension of the ascending node Ω(t) caused by orbital perturbations and a time-varying argument of perigee ω(t) as well as a change υ(t) in the true anomaly caused by the Earth's gravity; at each moment, in the Cartesian coordinate system with the geocenter as the origin, the position coordinates (x_sat, y_sat, z_sat) of the moving satellite are expressed as:

$x_{\mathrm{sat}}=R\left(t\right)·\left\{\cos \left(\omega \left(t\right)+v\left(t\right)\right)·\cos \Omega \left(t\right)-\sin \left(\omega \left(t\right)+v\left(t\right)\right)·\sin \Omega \left(t\right)·\cos \left(\iota \right)\right\}$ $y_{\mathrm{sat}}=R\left(t\right)·\left\{\cos \left(\omega \left(t\right)+v\left(t\right)\right)·\sin \Omega \left(t\right)-\sin \left(\omega \left(t\right)+v\left(t\right)\right)·\cos \Omega \left(t\right)·\cos \left(\iota \right)\right\}$ $z_{\mathrm{sat}}=R\left(t\right)·\sin \left(\omega \left(t\right)+v\left(t\right)\right)·\sin \left(\iota \right)$

• • where, R(t) is a distance from the geocenter to the satellite at each moment, calculated as:

$R\left(t\right)=\frac{a·{\left(1-e\right)}^2}{1+e·\cos v\left(t\right)}$

• • the in-orbit speed of a satellite in an elliptical orbit is time-varying, given by

$\sqrt{{\mu }_E\left(\frac{2}{R\left(t\right)}-\frac{1}{a}\right)}$

km/s; when the eccentricity e of the elliptical orbit is set to 0, the elliptical orbit is reduced to a circular orbit with a constant speed √{square root over (μ_E/a)} km/s, where, μ_E is a gravitational constant, with a value of 3.986012×10⁵ km³/s²; and

• • Step S202: the receiving end is set as a moving vehicle, of which a speed v_rx and a moving trajectory are set.

Further, the step S3 includes the following specific steps:

• • Step S301: generating spatially consistent large-scale parameters according to the frequency, environment and satellite elevation: large-scale fading PL, shadow fading SH, delay spread DS, azimuth angular spread of arrival ASA, elevation angular spread of arrival ESA, Rician K-factor KF, and cross-polarization ratio XPR, wherein for each large-scale parameter to be calculated, the following general formula is employed:

$V=V_{\mu }+V_{ϵ}·{\log }_{10}d+V_{\gamma }·{\log }_{10}{f}_{\mathrm{GHz}}+V_{\alpha }·{\log }_{10}{\alpha }_{\mathrm{rad}}+X\left(V_{\sigma }+V_{\delta }·{\log }_{10}{f}_{\mathrm{GHz}}+V_{\beta }·{\log }_{10}{\alpha }_{\mathrm{rad}}\right)$

• • where, the parameters V_μ, V_∈, V_γ, V_α, V_σ, V_δ and V_β are related to the environment, d is a distance between the transmitting end and the receiving end, f_GHz is a carrier frequency, α_rad is a satellite elevation, and X is a spatially consistent normally distributed random variable with a mean of 0 and a variance of 1; and
  • Step S302: setting a rainfall rate to R, wherein the effect of rainfall on channel multipath fading is reflected in the changes of the parameters KF, DS, ASA and ESA caused by the rainfall rate; the model employs the parameter ξ_KF to denote the effect of the rainfall rate on the Rician K-factor, and the Rician K-factor KF_R affected by the rainfall is expressed as:

${\mathrm{KF}}_R=\mathrm{KF}-R·{\xi }_{\mathrm{KF}}$

• • the model assumes that the effect of rainfall on multipath is linear, using the parameters ξ_DS, ξ_ASA and ξ_ESA to measure the effect of rainfall on a cluster; the delay spread DS_R and angular spread ASA_R and ESA_R affected by rainfall are expressed as:

${\mathrm{DS}}_R=\mathrm{DS}\left(1+R·{\xi }_{\mathrm{DS}}\right)$ ${\mathrm{ASA}}_R=\mathrm{ASA}\left(1+R·{\xi }_{\mathrm{ASA}}\right)$ ${\mathrm{ESR}}_R=\mathrm{ESA}\left(1+R·{\xi }_{\mathrm{ESA}}\right)$

• • in addition, an increase in the number of multipaths during rainfall is modeled as an increase in the number of clusters; and the model models the newly added clusters N_rain due to rainfall in a Poisson distribution:

$N_{\mathrm{rain}}\sim P\left(R·{\xi }_{\lambda }\right)$

• • where, P represents a Poisson distribution, and ξ_λ represents an expectation of the distribution.

Further, the step S4 includes the following specific steps:

• • Step S401: calculating a free path loss PL, and modeling as a log-distance path loss model, wherein the shadow fading SF follows a log-normal distribution;
  • Step S402: calculating the atmospheric absorption A_R, setting the receiving-end height equal to the sea level height, and taking the mean annual global reference atmosphere values as the environmental parameters: temperature T, dry air atmospheric pressure p, water vapor density ρ and water vapor pressure e, wherein A_G is expressed as:

$A_G=\frac{A_{\mathrm{zeni}\mathrm{th}}\left(f\right)}{\sin \left(\theta \right)}$

• • where, θ is the satellite elevation, and A_zenith(f) is the zenith angle attenuation; and
  • Step S403: calculating the large-scale fading caused by rainfall at a specific rainfall rate, and calculating a rainfall attenuation coefficient with reference to the Crane model, the rainfall attenuation expressed as:

$A_R=a_{\mathrm{Crane}}{R}^{b_{\mathrm{Crane}}}·\frac{h_R}{\sin {\theta }_{\mathrm{sat}}}$

• • where, h_R denotes the rainfall height, θ_sat denotes the satellite elevation, the parameters a_Crane and b_Crane are Crane model parameters which are obtained from discrete calculations and can be derived by curve fitting to a power law coefficient; and a relationship between the rainfall height h_R and the zero-degree isotherm h₀ is given by:

$h_R=h_0+0.36$

• • where, the unit thereof is km, and h₀ is related to the longitude and latitude of the earth.

Further, the step S5 specifically includes:

• • Step S501: calculating a delay τ_n of each cluster, a receiving-end azimuth angle φ_n^a and a receiving-end elevation angle θ_n^a using the large-scale parameters calculated in Step S301, wherein the delay of the cluster follows a unilateral exponential distribution, and an initial value thereof is calculated as:

${\tilde{\tau }}_l=-\ln \left\{X_l^{\tau }\left(x_t,y_t,z_t,x_r,y_r,z_r\right)\right\}$

• • where, (x_t, y_t, z_t) denotes the transmitting-end coordinates, and (x_r, y_r, z_r) denotes the receiving-end coordinates;
  • and scaled by the large-scale parameters as:

${\tau }_l={\tilde{\tau }}_l·\frac{1}{F}·\sum _{f=1}^F\frac{{\mathrm{DS}}_f}{\stackrel{\sim }{{\mathrm{DS}}_f}}$

• • where, f=1 . . . F denotes F carrier frequencies, _f is an initial delay spread, and DS_f is a large-scale parameter;
  • the initial value of angle ϕ_l follows a uniform distribution and is scaled by a large-scale parameter:

${\varphi }_l=\arg \exp \left\{j·{\tilde{\varphi }}_l·s\right\},s=\frac{1}{F}·\sum _{f=1}^F\frac{A{S}_f}{\stackrel{\sim }{{\mathrm{AS}}_f}},s<\{\begin{array}{c}3.,\mathrm{azimuth}\mathrm{angle}\\ 1.5,\mathrm{elevation}\mathrm{angle}\end{array}$

• • where, _f is an initial delay spread, AS_f is the large-scale parameter, and s is a scale;
  • Step S502: calculating the polar coordinates of each cluster with the receiving end as the origin via a geometrical relationship, according to the basic information of the cluster generated in Step S501:

${\left(d_n^T\right)}^2={\left(d_n^R\right)}^2+d^2-2·d_n^R·d·\cos {\alpha }_n$

• • where d_n^T is the distance from the nth cluster to the transmitting end at the initial moment, d_n^R is the distance from the nth cluster to the receiving end, and d is the distance between the transmitting end and the receiving end; the angle α_n is the angle between a unit vector of the receiving end pointing to the nth cluster and a unit vector of the receiving end pointing to the transmitting end; with the given elevation angle φ_n^a and θ_n^a of the nth scatterer and the three-dimensional angular information φ_sat and θ_sat of the transmitting end of the satellite, a unit vector ({circumflex over (x)}_n, ŷ_n, {circumflex over (z)}_n) of the receiving end pointing to the scatterer and a unit vector ({circumflex over (x)}_sat, ŷ_sat, {circumflex over (z)}_sat) of the receiving end pointing to the transmitting end can be calculated, and the angle α_n can be calculated by the following equation:

$\{\begin{array}{c}{\hat{x}}_n=\cos {\theta }_n^a\cos {\phi }_n^a\\ {\hat{y}}_n=\cos {\theta }_n^a\sin {\phi }_n^a\\ {\hat{z}}_n=\sin {\theta }_n^a\end{array},\{\begin{array}{c}{\hat{x}}_{\mathrm{sat}}=\cos {\theta }_{\mathrm{sat}}\cos {\phi }_{\mathrm{sat}}\\ {\hat{y}}_{\mathrm{sat}}=\cos {\theta }_{\mathrm{sat}}\sin {\phi }_{\mathrm{sat}}\\ {\hat{z}}_{\mathrm{sat}}=\sin {\theta }_{\mathrm{sat}}\end{array}$ $\cos {\alpha }_n=\frac{{\hat{x}}_n{\hat{x}}_{\mathrm{sat}}+{\hat{y}}_n{\hat{y}}_{\mathrm{sat}}+{\hat{z}}_n{\hat{z}}_{\mathrm{sat}}}{❘\left({\hat{x}}_n,{\hat{y}}_n,{\hat{z}}_n\right)❘❘\left({\hat{x}}_{\mathrm{sat}},{\hat{y}}_{\mathrm{sat}},{\hat{z}}_{\mathrm{sat}}\right)❘}$

• • and substituting α_n in the first equation in Step S502, to obtain the central position (d_n^R, φ_n^a, θ_n^a) of the cluster in the polar coordinates with the receiving end as the origin;
  • Step S503: calculating the geometric positions of the scatterers in each cluster, wherein the scatterers in the model are in a three-dimensional elliptical Gaussian distribution around the center of the cluster; the three dimensions of the scatterer follow a Gaussian distribution with a standard deviation of σ_x, σ_y, σ_z, and in the three-dimensional Cartesian coordinate system with the receiving end as the origin, the coordinates (x_nm, y_nm, z_nm) of each scatterer are calculated by the following equation:

$\left[\begin{array}{c}{x}_{n_m}\\ y_{n_m}\\ z_{n_m}\end{array}\right]=\left[\begin{array}{ccc}\cos \left({\phi }_n\right)& -\sin \left({\phi }_n\right)& 0\\ \sin \left({\phi }_n\right)& \cos \left({\phi }_n\right)& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{ccc}\cos \left({\theta }_n\right)& 0& -\sin \left({\theta }_n\right)\\ 0& 1& 0\\ \sin \left({\theta }_n\right)& 0& \cos \left({\theta }_n\right)\end{array}\right]·\left[\begin{array}{c}{x}^{\prime }+d_n^R\\ y^{\prime }\\ z^{\prime }\end{array}\right]$ $p\left(x^{\prime },y^{\prime },z^{\prime }\right)=\frac{\exp \left(-\frac{x^{\mathrm{\prime 2}}}{2{\sigma }_x^2}-\frac{y^{\mathrm{\prime 2}}}{2{\sigma }_y^2}-\frac{z^{\mathrm{\prime 2}}}{2{\sigma }_z^2}\right)}{{\left(2\pi \right)}^{3/2}{\sigma }_x{\sigma }_y{\sigma }_z}$

• • Step S504: calculating the path length d_qp,nm(t) and delay τ_qp,nm(t) of each path according to the geometric positions of the scatterers determined in Step S503, wherein the model needs to consider the delay resolution of subpaths, which is, with the changes of power in time domain and frequency domain taken into account, calculated as:

$P_{\mathrm{qp},n_m,f_c}\left(t,f\right)=\exp \left(-{\tau }_{\mathrm{qp},n_m}^{\prime }\left(t\right)g_{f_c}^{\mathrm{DS}}\right)10^{-\frac{Z_n}{10}}{\left(\frac{f}{f_c}\right)}^{\gamma }$

• • where, τ′_qp,nm (t) is a relative delay of each subpath, Z_n is a cluster shadow calculated for each cluster, which follows a Gaussian distribution with a mean of 0;

${\left(\frac{f}{f_c}\right)}^{\gamma }$

models the frequency-dependent characteristic of power in a millimeter-wave high-bandwidth channel, and γ is related to frequency; g_fc^DS is a delay spread scale factor, which is, with different calculation methods employed in single-frequency point modeling and multi-frequency point modeling, expressed as:

$g_{f_c}^{\mathrm{DS}}=\{\begin{array}{c}-1.5·\ln \left\{1.2·\frac{{\mathrm{DS}}_{f_c}}{\max \left({\mathrm{DS}}_{f_c}\right)+\min \left({\mathrm{DS}}_{f_c}\right)}-0.15\right\},\mathrm{multiple}\mathrm{frequency}\mathrm{points}\\ \frac{r_{\tau }-1}{r_{\tau }\mathrm{DS}},\mathrm{single}\mathrm{frequency}\mathrm{point}\end{array}$

• • where, DS_fc denotes the large-scale parameters at different frequencies, DS is the delay spread, and r_τ is a delay distribution scale factor; and
  • Step S505: calculating the channel coefficient according to the generated parameters:

$H={\left[\mathrm{PL}·\mathrm{SH}·A_G·A_R\right]}^{1/2}·h_{\mathrm{qp},f_c}\left(t,\tau \right)$

• • where, the large-scale fading [PL·SH·A_G·A_R]^1/2 is calculated in Step 4, and the small-scale fading h_qp,fc(t, τ) is expressed as:

$h_{\mathrm{qp},f_c}\left(t,\tau \right)=\sqrt{\frac{K_R\left(t\right)}{K_R\left(t\right)+1}{h}_{\mathrm{qp},f_c}^L\left(t,\tau \right)}+\sqrt{\frac{1}{K_R\left(t\right)+1}{h}_{\mathrm{qp},f_c}^N\left(t,\tau \right)}$

• • where, K_R(t) is the Rician K-factor that varies with time, h_qp,fc^L(t,τ) is the channel impulse response of a LOS path, and h_qp,fc^N(t,τ) is the channel impulse response of a NLOS path, respectively expressed as:

$h_{\mathrm{qp},f_c}^L\left(t,\tau \right)={\left[\begin{array}{c}{F}_{q,V}\left({\varphi }_{E,L}^R\left(t\right),{\varphi }_{A,L}^R\left(t\right)\right)\\ F_{q,H}\left({\varphi }_{E,L}^R\left(t\right),{\varphi }_{A,L}^R\left(t\right)\right)\end{array}\right]}^T·\left[\begin{array}{cc}{e}^{j{\theta }_L^{\mathrm{VV}}}& 0\\ 0& -e^{j{\theta }_L^{\mathrm{HH}}}\end{array}\right]·F_r·\left[\begin{array}{c}{F}_{p,V}\left({\varphi }_{E,L}^T\left(t\right),{\varphi }_{A,L}^T\left(t\right)\right)\\ F_{p,H}\left({\varphi }_{E,L}^T\left(t\right),{\varphi }_{A,L}^T\left(t\right)\right)\end{array}\right]·e^{j2\pi f_c{\tau }_{\mathrm{qp}}^L\left(t\right)}\delta \left(\tau -{\tau }_{\mathrm{qp}}^L\left(t\right)\right)$ $h_{\mathrm{qp},f_c}^N\left(t,\tau \right)=\sum _{n=1}^{N_{\mathrm{qp}}\left(t\right)}\sum _{m=1}^{M_n\left(t\right)}{\left[\begin{array}{c}{F}_{q,f_c,V}\left({\varphi }_{E,n_m}^R\left(t\right),{\varphi }_{A,n_m}^R\left(t\right)\right)\\ F_{q,f_c,H}\left({\varphi }_{E,n_m}^R\left(t\right),{\varphi }_{A,n_m}^R\left(t\right)\right)\end{array}\right]}^T·\left[\begin{array}{cc}{e}^{j{\theta }_{n_m}^{\mathrm{VV}}}& \sqrt{{\mathrm{\mu \kappa }}_{n_m}^{-1}\left(t\right)}{e}^{j{\theta }_{n_m}^{\mathrm{HV}}}\\ \sqrt{{\kappa }_{n_m}^{-1}\left(t\right)}{e}^{j{\theta }_{n_m}^{\mathrm{HV}}}& \sqrt{\mu }{e}^{j{\theta }_{n_m}^{\mathrm{HH}}}\end{array}\right]·F_r·\left[\begin{array}{c}{F}_{p,f_c,V}\left({\varphi }_{E,n_m}^T\left(t\right),{\varphi }_{A,n_m}^T\left(t\right)\right)\\ F_{p,f_c,H}\left({\varphi }_{E,n_m}^T\left(t\right),{\varphi }_{A,n_m}^T\left(t\right)\right)\end{array}\right]·\sqrt{P_{\mathrm{qp},n_m,f_c}\left(t\right)}·e^{j2\pi f_c{\tau }_{\mathrm{qp}},n_m\left(t\right)}\delta \left(\tau -{\tau }_{\mathrm{qp},n_m}\left(t\right)\right)$

• • where, [·]^T denotes transpose, F_p(q),fc^,V is a horizontal polarization antenna pattern of the transmitting and receiving antennas, and F_p(q),fc^,H is an horizontal polarization antenna pattern; κ_nm(t) is a cross-polarization ratio, μ denotes co-polarization imbalance, ϕ_A,L^T(t) is an azimuth angle of departure corresponding to the LOS path at the moment t stretching from A₁^T to A₁^R, ϕ_E,L^T(t) is an elevation angle of departure corresponding to the LOS path at the moment t stretching from A₁^T to A₁^R, and ϕ_A,L^R(t) is an azimuth angle of arrival corresponding to the LoS path at the moment t stretching from A₁^T to A₁^R; ϕ_E,L^R(t) is an elevation angle of arrival corresponding to the LoS path at the moment t stretching from A₁^T to A₁^R, θ_L^VV and θ_L^HH represent the initial phases of the LOS path, θ_nm^HH, θ_nm^HV, θ_nm^VV and θ_nm^VH represent the initial phases of the NLOS path, which are random variables that follow a uniform distribution between 0 and 2π; Faraday rotation matrix F_r refers to the rotation of a polarization plane caused by the propagation of electromagnetic waves through the ionosphere in the satellite scenario, which needs to be considered in communication scenarios below 10 GHZ; P_qp,nm,fc(t) is the subpath power calculated in Step S504, τ_qp^L(t) is the absolute delay of the Mth of LOS subpath of LOS, which can be calculated by dividing the distance between the transmitting and receiving antennas q,p by the speed of light, and τ_qp,nm(t) is the absolute delay of the mth subpath in the nth cluster, which can be calculated by dividing the length of the n_mth subpath between the transmitting and receiving antennas q, p by the speed of light.

Further, the step S6 includes the following specific steps:

• • Step S601: introducing a survival probability P_surv(Δt, Δf) of the cluster in the time-frequency domain, with the birth and death of the cluster in the time domain and frequency domain taken into account:

$P_{\mathrm{surv}}\left(\Delta t,\Delta f\right)=P_{\mathrm{surv}}\left(\mathrm{\Delta t}\right)·P_{\mathrm{surv}}\left(\Delta f\right)=e^{-{\lambda }_R\frac{v_{\mathrm{rx}}\Delta t}{D_c^T}}·e^{-{\lambda }_R\frac{\Delta f}{D_c^F}}$

• • where, P_surv(Δt) is the time-domain cluster survival probability, P_surv(Δf) is the frequency-domain cluster survival probability, Δt is a time interval, Δf is a frequency interval, the birth-death process of the cluster is described jointly by a generation rate λ_G and a disappearance rate λ_R of the cluster, these two parameters are related to the environmental characteristics of the communication scenario and the antenna pattern; parameters D_c^T, D_c^F are scenario-related factors in the time domain and frequency domain, obtained by channel measurement in a specific scenario; and an expected number {N_new} of new clusters is calculated as; and

$𝔼\left\{N_{\mathrm{new}}\right\}=\frac{{\lambda }_G}{{\lambda }_R}\left[1-P_{\mathrm{surv}}\left(\Delta t,\Delta f\right)\right]$

• • Step S602: updating the channel coefficient according to the step S5, with the movement of the transmitting and receiving ends and the clusters as well as the birth and death of the clusters taken into account.

Further, the step S7 specifically includes:

• • Step S701: calculating an effective path loss, which consists of three parts: free space path loss PL, atmospheric absorption A_G and rainfall attenuation A_R, wherein the three large-scale fading components do not change rapidly with time, showing a relatively constant trend, and the sum of the three is defined as the effective path loss: PL_eff=PL+A_G+A_R
  • Step S702: calculating the root-mean-square (RMS) delay spread of the channel, expressed as:

$D_{\mathrm{qp},\mathrm{rms}}\left(t\right)=\sqrt{\frac{{\int }_{-\infty }^{\infty }{\left(\tau -{\mu }_{\mathrm{qp},\tau }\left(t\right)\right)}^2{P}_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau }{{\int }_{-\infty }^{\infty }{P}_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau }}$

• • where, τ is the delay, and the average delay can be expressed

${\mu }_{\mathrm{qp},\tau }\left(t\right)=\frac{{\int }_{-\infty }^{\infty }\tau ·P_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau }{{\int }_{-\infty }^{\infty }{P}_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau };$

• • Step S703: calculating a time-frequency correlation function of the channel, with a theoretical value expressed as:

$R_{qp}\left(t,f;\Delta t,\Delta f\right)=𝔼\left[H_{qp}\left(t,f\right)H_{qp}^{*}\left(t+\Delta t,f+\Delta f\right)\right]$

• • where, [·] represents the expectation, [·]* represents the conjugate of a complex number, H_qp represents a channel transfer function, which is obtained by Fourier transform of the channel impulse response; without considering the correlation between the LOS path and the NLOS path, the time-frequency correlation function is expressed as the sum of the time-frequency correlation functions of the paths:

$R_{qp}\left(t,f;\Delta t,\Delta f\right)=R_{qp}^L\left(t,f;\Delta t,\Delta f\right)+R_{qp}^N\left(t,f;\Delta t,\Delta f\right)$

• • where, R_qp^L(t, f; Δt, Δf) represents the time-frequency correlation function of the LOSpath, expressed as:

$R_{qp}^L\left(t,f;\Delta t,\Delta f\right)=𝔼[\sqrt{\frac{K_R\left(t\right)}{K_R\left(t\right)+1}·\frac{K_R\left(t+\Delta t\right)}{K_R\left(t+\Delta t\right)+1}}{e}^{j2\pi (\frac{\left(f_c-f\right)d_{\mathrm{qp}}\left(t\right)-\left(f_c-f-\Delta f\right)\mathrm{dqp}}{c}\left(t+\Delta t\right)})]$

• • where, c is the speed of light;
  • R_qp^N(t, f; Δt, Δf) represents the time-frequency correlation function of the NLOS path, and without considering the correlation between multipaths, is expressed as:

$R_{qp}^N\left(t,f;\Delta t,\Delta f\right)=P_{\mathrm{surv}}\left(\Delta t,\Delta f\right)·𝔼[\sum _{m=1}^M\sum _{n=1}^N\sqrt{P_{\mathrm{qp},n_m,f_c}\left(t,f\right)P_{\mathrm{qp},n_m,f_c}\left(t+\Delta t,f+\Delta f\right)}·e^{j2\pi \left(\frac{\left(f_c-f\right)d_{\mathrm{qp},}{n}_m\left(t\right)-\left(f_c-f-\Delta f\right)d_{\mathrm{qp},}{n}_m\left(t+\Delta t\right)}{c}\right)}|$

• • Step S704: calculating the Doppler frequency of each multipath:

$v_{\mathrm{qp},n_m}\left(t\right)=-\frac{v^T}{\lambda }\cos \left({\omega }_p^T\right)-\frac{v^R}{\lambda }\cos \left({\omega }_q^R\right)+\frac{{\sin }^2\left({\omega }_p^T\right){\left(v^T\right)}^{2}t}{\lambda \left(d_{n_m}^T-\cos \left({\vartheta }^T\right){\delta }_p\right)}+\frac{{\sin }^2\left({\omega }_q^R\right){\left(v^R\right)}^{2}t}{\lambda \left(d_{n_m}^R-\cos \left({\vartheta }^R\right){\delta }_q\right)}$

• • where δ_p/δ_q in the equation represent the antenna spacings at the transmitting end and the receiving end, ω_p^T represents an angle between the moving direction of the transmitting end and the n_mth path corresponding to the pth transmitting antenna, ω_q^R represents an angle between the moving direction of the receiving end and the n_mth path corresponding to the qth receiving antenna, and ϑ^T represents an angle between the transmitting antenna array and the n_mth path corresponding to the 1st transmitting antenna/receiving antenna; and by calculating the Doppler frequencies corresponding to all multipaths in each simulation, local Doppler spread under this simulation can be obtained, and through multiple simulations to obtain the sample mean, local Doppler spread under this type of scenario can be obtained:

$B_{qp}^{\left(2\right)}\left(t\right)=𝔼\left[\sqrt{\stackrel{\_}{\left[{v_{\mathrm{qp},n_m}\left(t\right)}^2\right]}-\stackrel{\_}{{\left[v_{\mathrm{qp},n_m}\left(t\right)\right]}^2}}\right].$

The satellite communication-oriented geometry-based stochastic channel modeling method of the present invention has the following advantages: The present invention provides a general three-dimensional geometry-based stochastic satellite channel model, which is applicable to the commonly used satellite communication bands ranging from S-band to millimeter-wave band, and, by modeling the satellite mobility based on real orbits, is suitable for the communication scenarios of LEO, MEO and GEO satellites. On the basis of a satellite channel simulation scenario, the present invention sets the scenario layout parameters, initializes the trajectories and speed of the satellite and the receiving end, and supports simultaneous movement of the transmitting and receiving ends and the clusters. It takes into account the spatial consistency of large-scale parameters and the correlation of the channel parameters: environment, frequency and satellite elevation, and considers the effect of rainfall on the large-scale parameters. With respect to large-scale fading, it considers path loss, shadow fading, atmospheric absorption, and rainfall attenuation. The time delay, angle and power of the cluster are calculated by calculating the geometric position information of the cluster and the scatterer, to generate the channel coefficient. The birth and death of the cluster in the time domain and the frequency domain are taken into account, and the statistical characteristics of the channel are derived for simulation analysis. In addition, the model considers the ionospheric effect that is more likely to happen at the low-frequency band and the rainfall attenuation that is more likely to happen at the millimeter-wave band respectively, becoming the first geometry-based stochastic channel model that considers the effect of rainfall attenuation on both the large-scale fading and small-scale fading.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic flow chart of a satellite communication-oriented geometry-based stochastic channel modeling method according to the embodiment; and

FIG. 2 is a schematic diagram of a three-dimensional general geometry-based stochastic channel model for satellite communication according to the embodiment.

DETAILED DESCRIPTION OF THE INVENTION

In order to better understand the purpose, structure and function of the present invention, the satellite communication-oriented geometry-based stochastic channel modeling method of the present invention is further described in detail below with reference to the accompanying drawings.

Embodiment 1

Referring to FIG. 1 and FIG. 2, this embodiment provides a scenario prediction channel modeling method based on scatterer density, specifically including:

• • Step S1: a satellite channel simulation scenario is established, and the physical environment parameters in the satellite communication channel scenario are set, as shown in FIG. 2.

Specifically, in this embodiment, the step S1 specifically includes:

• • Step S101: the transmitting end and receiving end of the model are a satellite mobile receiving end and a ground mobile receiving end respectively, the orbital altitude h_sat of the satellite needs to be pre-set, and the environment of the ground mobile receiving end also needs to be pre-set as a dense city, a city, a suburb, or a rural area;
  • Step S102: both the transmitting end and receiving end of the model are a multiple-transmit multiple-receive antenna array, so the spacing and angle parameters of the antenna array need to be initialized, where the antenna array is a uniformly distributed linear array. The transmitting end has a total of P antennas, A_p^T denotes the pth antenna, and β_A^T and β_E^T denote an azimuth angle and an elevation angle of the transmitting-end antenna array respectively, a spacing between the antennas is δ^T; and the receiving end has a total of Q antennas, A_q^T denotes the qth antenna, β_A^R and β_E^R denote an azimuth angle and an elevation angle of the receiving-end antenna array respectively, and a spacing between the antennas is δ^R.
  • Step S103: the system parameters of the model further include a carrier frequency f_c, and a rainfall rate set based on whether it is raining, where, without rain, it is set to 0 mm/h; and with a heavy rain, it is set to 50 mm/h.
  • Step S2: initializing the trajectory and speed of a satellite and a receiving end;

Specifically, in this embodiment, the step S2 specifically includes:

• • Step S201: the initial position of the satellite can be determined together by the orbital altitude h_sat, elevation θ_sat and azimuth φ_sat. Satellites can be divided into GEO satellites and NGEO satellites depending on whether the satellites are stationary relative to the earth. GEO satellites are stationary relative to the ground at an altitude of about 36,000 km. Due to the Earth's gravity, NGEO satellites orbit the Earth along a certain trajectory. First, an elliptical orbital plane of the satellite is determined by setting the following five parameters: a length a of the semi-major axis of the elliptical orbit; an eccentricity e of the elliptical orbit, where, if the eccentricity e is set to 0, the elliptical orbit is reduced to a circular orbit; an angle between the orbit and the equatorial plane is an orbital plane inclination ι; the longitude of the point where the orbit crosses the equatorial plane upward is called a right ascension of the ascending node Ω; the direction of the ellipse on the orbital plane is determined by an argument of perigee ω of an angle between the orbital perigee and the ascending node. A true anomaly υ represents an angle swept by moving along the orbit from the perigee; and after determination of the orbital plane, the specific position of the satellite in the elliptical orbit can be determined via υ; The satellite trajectory is determined together by a time-varying right ascension of the ascending node Ω(t) caused by orbital perturbations and a time-varying argument of perigee ω(t) as well as a change υ(t) in the true anomaly caused by the Earth's gravity, where the three values are calculated by referring to a standardization document QuaDRIGA. At each moment, in the Cartesian coordinate system with the geocenter as the origin, the position of the moving satellite can be expressed as:

$x_{sat}=R\left(t\right)·\left\{\cos \left(\omega \left(t\right)+v\left(t\right)\right)·\cos \Omega \left(t\right)-\sin \left(\omega \left(t\right)+v\left(t\right)\right)·\sin \Omega \left(t\right)·\cos \left(\iota \right)\right\}$ $y_{sat}=R\left(t\right)·\left\{\cos \left(\omega \left(t\right)+v\left(t\right)\right)·\sin \Omega \left(t\right)-\sin \left(\omega \left(t\right)+v\left(t\right)\right)·\cos \Omega \left(t\right)·\cos \left(\iota \right)\right\}$ $z_{sat}=R\left(t\right)·\sin \left(\omega \left(t\right)+v\left(t\right)\right)·\sin \left(\iota \right)$

• • where, R(t) is a distance from the geocenter to the satellite at each moment, calculated as:

$R\left(t\right)=\frac{a·{\left(1-e\right)}^2}{1+e·\cos v\left(t\right)}$

• • the in-orbit speed of a satellite in an elliptical orbit is time-varying, given by

$\sqrt{{\mu }_E\left(\frac{2}{R\left(t\right)}-\frac{1}{a}\right)}$

km/s; when the eccentricity e of the elliptical orbit is set to 0, the elliptical orbit is reduced to a circular orbit with a constant speed √{square root over (μ_E/a)} km/s, where, μ_E is a gravitational constant, with a value of 3.986012×10⁵ km³/s².

• • Step S202: the receiving end is set as a moving vehicle, of which a speed v_rx and a moving trajectory are set.
  • Step S3: calculating spatially consistent large-scale parameters, and the effect of rainfall on the large-scale parameters;

Specifically, in this embodiment, step S3 specifically includes:

• • Step S301: in order to generate the large-scale parameters related to frequency, environment and elevation, this model adopts the linear model and satellite scenario parameter table proposed in the standardization document QuaDRIGA and related studies thereof to generate spatially consistent large-scale parameters: large-scale fading PL, shadow fading SH, delay spread DS, azimuth angular spread of arrival ASA, elevation angular spread of arrival ESA, Rician K-factor KF and cross-polarization ratio XPR. For each large-scale parameter that needs to be calculated, the following general formula is employed:

$V=V_{\mu }+V_{ϵ}·{\log }_{10}d+V_{\gamma }·{\log }_{10}{f}_{GHz}+V_a·{\log }_{10}{\alpha }_{\mathrm{rad}}+X\left(V_{\sigma }+V_{\delta }·{\log }_{10}{f}_{GHz}+V_{\beta }·{\log }_{10}{\alpha }_{\mathrm{rad}}\right)$

• • where, the parameters such as V_μ are related to the environment, and X is a spatially consistent normally distributed random variable with a mean of 0 and a variance of 1.

Step S302: a rainfall rate R is considered, where the effect of rainfall on channel multipath fading is reflected in the changes of the parameters KF, DS, ASA and ESA caused by the rainfall rate. The existing GBSM satellite channel model assumes that the tropospheric effect and multipath effect are independent of each other. In fact, rainfall may cause atmospheric inhomogeneity and consequent strong and compact edges of rain cells, which in turn lead to stronger multipath effects. In addition, rain adheres to the surfaces of scatterers such as buildings and trees, causing the scatterers to change the properties such as scattering and reflection of electromagnetic waves. Through actual measurements, it is found that rainfall will increase the number of multipaths, and as the rainfall rate rises, the delay spread and power of the multipaths increase, while the power of the direct path decreases. The Rician K-factor KF decreases as the rainfall rate rises, and the multipath effect of the channel becomes more obvious. This model employs the parameter ξ_KF to denote the effect of the rainfall rate on the Rician K-factor, and the Rician K-factor affected by the rainfall can be modeled as:

$K{F}_R=KF-R·{\xi }_{KF}$

As the rainfall rate increases, the multipath delay spread increases; and due to the presence of compact rain cells and the change in scatterer humidity, new multipaths appear, and it can be inferred that rainfall will cause a certain increase in angular spread. This model assumes that the effect of rainfall on multipath is linear, using the parameters ξ_DS, ξ_ASA and ξ_ESA to measure the effect of rainfall on a cluster; and the specific value needs to be determined through measurement. The effect of the rainfall rate on the delay spread and angular spread is expressed as:

${\mathrm{DS}}_R=\mathrm{DS}\left(1+R·{\xi }_{\mathrm{DS}}\right)$ ${\mathrm{ASA}}_R=\mathrm{ASA}\left(1+R·{\xi }_{\mathrm{ASA}}\right)$ ${\mathrm{ESA}}_R=\mathrm{ESA}\left(1+R·{\xi }_{\mathrm{ESA}}\right)$

• • in addition, an increase in the number of multipaths during rainfall may be modeled as an increase in the number of clusters; and the model models the newly added clusters N_rain due to rainfall in a Poisson distribution:

$N_{\mathrm{rain}}\sim P\left(R·{\xi }_{\lambda }\right)$

• • Step S4: calculating a path loss, a shadow fading, an atmospheric absorption and a rainfall attenuation;

Specifically, in this embodiment, the step S4 specifically includes:

• • Step S401: calculating a free path loss PL, and modeling as a log-distance path loss model, where the shadow fading SF follows a log-normal distribution;
  • Step S402: calculating the atmospheric absorption A_R, by referring to the calculation method of approximate estimation for gas attenuation given in the Recommendations ITU-R-P.676, assuming that the receiving-end height is equal to the sea level height, and taking the mean annual global reference atmosphere values as the environmental parameters: temperature T=288.15 K, dry air atmospheric pressure p=1013.25 hPa, water vapor density ρ=7.5 g/m³ and water vapor pressure

$e=\frac{\rho T}{216.5}=9.98\mathrm{hPa}.$

A_G may be expressed as:

$A_G=\frac{A_{\mathrm{zenith}}\left(f\right)}{\sin \left(\theta \right)}$

• • where, θ is the satellite elevation, and A_zenith(f) is the zenith angle attenuation, depending on the environmental parameters and carrier frequency; and
  • Step S403: calculating the large-scale fading caused by rainfall at a specific rainfall rate, and calculating a rainfall attenuation coefficient with reference to the Crane model, the rainfall attenuation expressed as:

$A_R=a_{\mathrm{Crane}}{R}^{b_{\mathrm{Crane}}}·\frac{h_R}{\sin {\theta }_{\mathrm{sat}}}$

• • where, h_R denotes the rainfall height, θ_sat denotes the satellite elevation, the parameters a_Crane and b_Crane are Crane model parameters which are obtained from discrete calculations and can be derived by curve fitting to a power law coefficient; and a relationship between the rainfall height h_R and the zero-degree isotherm h₀ is given by:

$h_R=h_0+0.36$

• • where, the unit thereof is km, and h₀ is related to the longitude and latitude of the earth. where, the parameters a_Crane and b_Crane are calculated as:

${\log }_{10}{a}_{\mathrm{Crane}}=\sum _{j=1}^4\left({\alpha }_j\exp \left[-{\left(\frac{{\log }_{10}{f}_c-{\beta }_j}{{\gamma }_j}\right)}^2\right]\right)+m_k{\log }_{10}{f}_c+c_k{b}_{\mathrm{Crane}}=\sum _{j=1}^5({\alpha }_j\exp \left[-{\left(\frac{{\log }_{10}{f}_c-{\beta }_j}{{\gamma }_j}\right)}^2\right]+m_{\alpha }{\log }_{10}{f}_c+c_{\alpha }$

• • where, the values of the parameters α, β, γ, m are obtained by referring to Table 1-4 in Recommendations ITU-R P.838-3. Note: h_R is not the true altitude of the satellite because rainfall is only distributed within a certain altitude in the troposphere, referring to the relationship between the rainfall height h_R and the zero-degree isotherm h₀ given in Recommendations ITU-R P.839-4:

$h_R=h_0+0.36$

• • the unit thereof is km, and h₀ is related to the longitude and latitude of the earth. Since this model does not consider the spatial statistical characteristics of the Earth's scale, so take the mean value to be 3 km. If necessary, accurate h₀ may also be obtained according to the digital map given in the above standardization documents.
  • Step S5: initializing the central positions of a cluster and a scatterer, and calculating the delay, angle and power of the cluster according to the geometric position information of the transmitting and receiving ends and the scatterer, to generate a channel coefficient;

Specifically, in this embodiment, the step S5 specifically includes:

• • Step S501: calculating a delay τ_n of each cluster, a receiving-end azimuth angle φ_n^a and a receiving-end elevation angle θ_n^a using the large-scale parameters calculated in Step S301, wherein the delay of the cluster follows a unilateral exponential distribution, and an initial value thereof is calculated as:

${\tilde{\tau }}_l=-\ln \left\{X_l^{\tau }\left(x_t,y_t,z_t,x_r,y_r,z_r\right)\right\}$

• • and may be scaled by the large-scale parameters as:

${\tau }_l={\tilde{\tau }}_l·\frac{1}{F}·\sum _{f=1}^F\frac{D{S}_f}{}$

• • the initial value of angle {tilde over (ϕ)}_l follows a uniform distribution and may be scaled by a large-scale parameter:

${\varphi }_l=\arg \exp \left\{j·{\tilde{\varphi }}_l·s\right\},s=\frac{1}{F}·\sum _{f=1}^F\frac{A{S}_f}{},s<\{\begin{array}{cc}3.,& \mathrm{azimuth}\mathrm{angle}\\ 1.5,& \mathrm{elevation}\mathrm{angle}\end{array}$

• • Step S502: calculating the polar coordinates of each cluster with the receiving end as the origin via a geometrical relationship, according to the basic information of the cluster generated in Step S501:

${\left(d_n^T\right)}^2={\left(d_n^R\right)}^2+d^2-2·d_n^R·d·\cos {\alpha }_n$

• • where, d_n^T is the distance from the nth cluster to the transmitting end at the initial moment, d_n^R is the distance from the nth cluster to the receiving end, and d is the distance between the transmitting end and the receiving end. The angle α_n is the angle between a unit vector of the receiving end pointing to the nth cluster and a unit vector of the receiving end pointing to the transmitting end; with the given three-dimensional angular information φ_n^a, θ_n^a, φ_sat, θ_sat of the nth scatterer and the transmitting end of satellite, the angle α_n may be calculated by the following equation:

$$\begin{gathered} \{\begin{array}{c}{\hat{x}}_n=\cos {\theta }_n^a\cos {\phi }_n^a\\ {\hat{y}}_n=\cos {\theta }_n^a\sin {\phi }_n^a\\ {\hat{z}}_n=\sin {\theta }_n^a\end{array},\{\begin{array}{c}{\hat{x}}_{\mathrm{sat}}=\cos {\theta }_{\mathrm{sat}}\cos {\phi }_{\mathrm{sat}}\\ {\hat{y}}_{\mathrm{sat}}=\cos {\theta }_{\mathrm{sat}}\sin {\phi }_{\mathrm{sat}}\\ {\hat{z}}_{\mathrm{sat}}=\sin {\theta }_{\mathrm{sat}}\end{array},\text{ \\ cos⁢αn=xˆn⁢xˆsat+yˆn⁢yˆsat+zˆn⁢zˆsat❘"[LeftBracketingBar]"(xˆn,yˆn,zˆn)❘"[RightBracketingBar]"⁢❘"[LeftBracketingBar]"(xˆsat,yˆsat,zˆsat)❘"[RightBracketingBar]"} \end{gathered}$$

• • next, the central position (d_n^R, φ_n^a, θ_n^a) of the cluster in the polar coordinates with the receiving end as the origin may be obtained;
  • Step S503: calculating the geometric positions of the scatterers in each cluster, where the scatterers in the model are in a three-dimensional elliptical Gaussian distribution around the center of the cluster; the three dimensions of the scatterer follow a Gaussian distribution with a standard deviation of σ_x, σ_y, σ_z, and in the three-dimensional Cartesian coordinate system with the receiving end as the origin, the coordinates (x_nm, y_nm, z_nm) of each scatterer may be calculated by the following equation:

$\left[\begin{array}{c}{x}_{n_m}\\ y_{n_m}\\ z_{n_m}\end{array}\right]=\left[\begin{array}{ccc}\cos \left({\phi }_n\right)& -\sin \left({\phi }_n\right)& 0\\ \sin \left({\phi }_n\right)& \cos \left({\phi }_n\right)& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{ccc}\cos \left({\theta }_n\right)& 0& -\sin \left({\theta }_n\right)\\ 0& 1& 0\\ \sin \left({\theta }_n\right)& 0& \cos \left({\theta }_n\right)\end{array}\right]·\left[\begin{array}{c}{x}^{\prime }+d_n^R\\ y^{\prime }\\ z^{\prime }\end{array}\right]p\left(x^{\prime },y^{\prime },z^{\prime }\right)=\frac{\exp \left(-\frac{x^{\mathrm{\prime 2}}}{2{\sigma }_x^2}-\frac{y^{\mathrm{\prime 2}}}{2{\sigma }_y^2}-\frac{z^{\mathrm{\prime 2}}}{2{\sigma }_z^2}\right)}{{\left(2\pi \right)}^{3/2}{\sigma }_x{\sigma }_y{\sigma }_z}$

• • Step S504: calculating the path length d_qp,nm(t), delay τ_qp,nm(t) and power of each path according to the geometric positions of the scatterers determined in Step S503, wherein the model needs to consider the delay resolution of subpaths, which, with the changes of power in the time domain and frequency domain taken into account, may be calculated as:

$P_{\mathrm{qp},n_m,f_c}\left(t,f\right)=\exp \left(-{\tau }_{\mathrm{qp},n_m}^{\prime }\left(t\right)g_{f_c}^{\mathrm{DS}}\right)10^{-\frac{Z_n}{10}}{\left(\frac{f}{f_c}\right)}^{\gamma }$

• • where, τ_qp,nm′(t) is a relative delay of each subpath, Z_n is a cluster shadow calculated for each cluster, which follows a Gaussian distribution with a mean of 0;

${\left(\frac{f}{f_c}\right)}^{\gamma }$

models the frequency-dependent characteristic of power in a millimeter-wave high-bandwidth channel, and γ is related to frequency; g_Fc^DS is a delay spread scale factor, which is, with different calculation methods employed in single-frequency point modeling and multi-frequency point modeling, expressed as:

$g_{f_c}^{\mathrm{DS}}=\{\begin{array}{c}-1.5·\ln \left\{1.2·\frac{{\mathrm{DS}}_{f_c}}{\max \left({\mathrm{DS}}_{f_c}\right)+\min \left({\mathrm{DS}}_{f_c}\right)}-0.15\right\},\mathrm{multiple}\mathrm{frequency}\mathrm{points}\\ \frac{r_{\tau }-1}{r_{\tau }\mathrm{DS}},\mathrm{single}\mathrm{frequency}\mathrm{point}\end{array}$

• • where, DS is the delay spread, and r_τ is a delay distribution scale factor; and
  • Step S505: calculating the channel coefficient according to the generated parameters:

$H={\left[\mathrm{PL}·\mathrm{SH}·A_G·A_R\right]}^{1/2}·h_{\mathrm{qp},f_c}\left(t,\tau \right)$

• • where, the large-scale fading [PL·SH·A_G·A_R]^1/2 is calculated in Step 4, and the small-scale fading h_qp,fc(t, τ) is expressed as:

$h_{\mathrm{qp},f_c}\left(t,\tau \right)=\sqrt{\frac{K_R\left(t\right)}{K_R\left(t\right)+1}{h}_{\mathrm{qp},f_c}^L\left(t,\tau \right)}+\sqrt{\frac{1}{K_R\left(t\right)}{h}_{\mathrm{qp},f_c}^N\left(t,\tau \right)}$

• • where, K_R(t) is the Rician K-factor that varies with time, and h_qp,fc^L(t, τ) and h_qp,fc^N(t, τ) are the channel impulse response of a LOS path and a NLOS path respectively, expressed as:

$\begin{array}{c}{h}_{\mathrm{qp},f_c}^L\left(t,\tau \right)={\left[\begin{array}{c}{F}_{q,V}\left({\varphi }_{E,L}^R\left(t\right),{\varphi }_{A,L}^R\left(t\right)\right)\\ F_{q,H}\left({\varphi }_{E,L}^R\left(t\right),{\varphi }_{A,L}^R\left(t\right)\right)\end{array}\right]}^T·\left[\begin{array}{cc}{e}^{j{\theta }_L^{\mathrm{VV}}}& 0\\ 0& -e^{j{\theta }_L^{\mathrm{HH}}}\end{array}\right]·F_r·\left[\begin{array}{c}{F}_{p,V}\left({\varphi }_{E,L}^T\left(t\right),{\varphi }_{A,L}^T\left(t\right)\right)\\ F_{p,H}\left({\varphi }_{E,L}^T\left(t\right),{\varphi }_{A,L}^T\left(t\right)\right)\end{array}\right]·e^{j2\pi f_c{\tau }_{\mathrm{qp}}^L\left(t\right)}\delta \left(\tau -{\tau }_{\mathrm{qp}}^L\left(t\right)\right)\\ h_{\mathrm{qp},f_c}^N\left(t,\tau \right)=\sum _{n=1}^{N_{\mathrm{qp}}\left(t\right)}\sum _{m=1}^{M_n\left(t\right)}{\left[\begin{array}{c}{F}_{q,f_c,,V}\left({\varphi }_{E,n_m}^R\left(t\right),{\varphi }_{A,n_m}^R\left(t\right)\right)\\ F_{q,f_c,H}\left({\varphi }_{E,n_m}^R\left(t\right),{\varphi }_{A,n_m}^R\left(t\right)\right)\end{array}\right]}^T·\left[\begin{array}{cc}{e}^{j{\theta }_{n_m}^{\mathrm{VV}}}& \sqrt{\mu {\kappa }_{n_m}^{-1}\left(t\right)}{e}^{j{\theta }_{n_m}^{\mathrm{VH}}}\\ \sqrt{{\kappa }_{n_m}^{-1}\left(t\right)}{e}^{j{\theta }_{n_m}^{\mathrm{HV}}}& \sqrt{\mu }{e}^{j{\theta }_{n_m}^{\mathrm{HH}}}\end{array}\right]·F_r·\left[\begin{array}{c}{F}_{p,f_c,,V}\left({\varphi }_{E,n_m}^T\left(t\right),{\varphi }_{A,n_m}^T\left(t\right)\right)\\ F_{p,f_c,H}\left({\varphi }_{E,n_m}^T\left(t\right),{\varphi }_{A,n_m}^T\left(t\right)\right)\end{array}\right]·\sqrt{P_{\mathrm{qp},n_m,f_c}\left(t\right)}·e^{j2\pi f_c{\tau }_{\mathrm{qp}},n_m\left(t\right)}\delta \left(\tau -{\tau }_{\mathrm{qp},n_m}\left(t\right)\right)\end{array}$

• • where, [.]^T denotes transpose, F_p(q),fcV is a vertical polarization antenna pattern of the transmitting and receiving antennas, and F_p(q),fcH is an horizontal polarization antenna pattern; κ_nm (t) is a cross-polarization ratio, u denotes co-polarization imbalance, ϕ_A,L^T(t) is an azimuth angle of departure corresponding to the LOS path at the moment t stretching from A₁^T to A₁^R, ϕ_E,L^T(t) is an elevation angle of departure corresponding to the LOS path at the moment t stretching from A₁^T to A₁^R, ϕ_A,L^R(t) is an azimuth angle of arrival corresponding to the LoS path at the moment t stretching from A₁^T to A₁^R; ϕ_E,L^R(t) is an elevation angle of arrival corresponding to the LoS path at the moment t stretching from A₁^T to A₁^R, θ_L^VV and θ_L^HH represent the initial phases of the LOS path, θ_nm^HH, θ_nm^HV, θ_nm^VV, and θ_nm^VH represent the initial phases of the NLOS path, which are random variables that follow a uniform distribution between 0 and 2π; Faraday rotation matrix F_r refers to the rotation of a polarization plane caused by the propagation of electromagnetic waves through the ionosphere in the satellite scenario, which needs to be considered in communication scenarios below 10 GHZ; P_qp,nm,fc(t) is the subpath power calculated in Step S504, τ_qp^L(t) is the absolute delay of the Mth of LOS subpath of LOS, which can be calculated by dividing the distance between the transmitting and receiving antennas q,p by the speed of light, and τ_qp,nm(t) is the absolute delay of the Mth subpath in the Nth cluster, which can be calculated by dividing the length of the n_mth subpath between the transmitting and receiving antennas q, p by the speed of light.

Step S6: updating the large-scale and small-scale parameters according to the movement of the transmitting and receiving ends and the birth-death process of the cluster, to generate a new channel coefficient; and

Specifically, in this embodiment, the step S6 specifically includes:

• • Step S601: introducing a survival probability P_surv(Δt, Δf) of the cluster in the time-frequency domain, with the birth and death of the cluster in the time domain and frequency domain taken into account:

$P_{\mathrm{surv}}\left(\Delta t,\Delta f\right)=P_{\mathrm{surv}}\left(\Delta t\right)·P_{\mathrm{surv}}\left(\Delta f\right)=e^{-{\lambda }_R\frac{v_{\mathrm{rx}}\Delta t}{D_c^T}}·e^{-{\lambda }_R\frac{\Delta f}{D_c^F}}$

• • where, P_surv(Δt) is the time-domain cluster survival probability, P_surv(Δf) is the frequency-domain cluster survival probability, Δt is a time interval, Δf is a frequency interval, the birth-death process of the cluster is described jointly by a generation rate λ_G and a disappearance rate λ_R of the cluster, these two parameters are related to the environmental characteristics of the communication scenario and the antenna pattern; parameters D_c^T, D_c^F are scenario-related factors in the time domain and frequency domain, obtained by channel measurement in a specific scenario; and an expected number {N_new} of new clusters is calculated as; and

$𝔼\left\{N_{\mathrm{new}}\right\}=\frac{{\lambda }_G}{{\lambda }_R}\left[1-P_{\mathrm{surv}}\left(\Delta t,\Delta f\right)\right]$

• • Step S602: updating the channel coefficient according to the step S5, with the movement of the transmitting and receiving ends and the clusters as well as the birth and death of the clusters taken into account.

Step S7: the statistical characteristics of the channel are derived, and simulation analysis is performed;

Specifically, in this embodiment, the step S7 specifically includes:

• • Step S701: calculating an effective path loss, which consists of three parts: free path loss PL, atmospheric absorption A_G and rainfall attenuation A_R, wherein the three large-scale fading components do not change rapidly with time, showing a relatively constant trend, and the sum of the three is defined as the effective path loss: PL_eff=PL+A_G+A_R
  • Step S702: calculating the root-mean-square (RMS) delay spread of the channel, expressed as:

$D_{\mathrm{qp},\mathrm{rms}}\left(t\right)=\sqrt{\frac{{\int }_{-\infty }^{\infty }{\left(\tau -{\mu }_{\mathrm{qp},\tau }\left(t\right)\right)}^2{P}_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau }{{\int }_{-\infty }^{\infty }{P}_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau }}$

• • where, τ is the delay, and the average delay can be expressed

${\mu }_{\mathrm{qp},\tau }\left(t\right)=\frac{{\int }_{-\infty }^{\infty }\tau ·P_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau }{{\int }_{-\infty }^{\infty }{P}_{\mathrm{qp},n_m,f_c}\left(t,\tau \right)d\tau };$

• • Step S703: calculating a time-frequency correlation function of the channel, with a theoretical value expressed as:

$R_{\mathrm{qp}}\left(t,f;\Delta t,\Delta f\right)=𝔼\left[H_{\mathrm{qp}}\left(t,f\right)H_{qp}^{*}\left(t+\Delta t,f+\Delta f\right)\right]$

• • where, [·] represents expectation, [·]* represents the conjugate of a complex number, H_qp represents a channel transfer function, which is obtained by Fourier transform of the channel impulse response; without considering the correlation between the LOS path and the NLOS path, the time-frequency correlation function is expressed as the sum of the time-frequency correlation functions of the paths:

$R_{\mathrm{qp}}\left(t,f;\Delta t,\Delta f\right)=R_{qp}^L\left(t,f;\Delta t,\Delta f\right)+R_{qp}^N\left(t,f;\Delta t,\Delta f\right)$

• • where, R_qp^L(t, f; Δt, Δf) represents the time-frequency correlation function of the LOS path, expressed as:

$R_{qp}^L\left(t,f;\Delta t,\Delta f\right)=𝔼\left[\sqrt{\frac{K_R\left(t\right)}{K_R\left(t\right)+1}·\frac{K_R\left(t+\Delta t\right)}{K_R\left(t+\Delta t\right)+1}{e}^{j2\pi \left(\frac{\left(f_c-f\right)d_{\mathrm{qp}}\left(t\right)-\left(f_c-f-\Delta f\right)\mathrm{dqp}}{c}\left(t+\Delta t\right)\right)}}\right]$

• • where, c is the speed of light;
  • R_qp^N(t, f; Δt, Δf) represents the time-frequency correlation function of the NLOS path, and without considering the correlation between multipaths, is expressed as:

$R_{qp}^N\left(t,f;\Delta t,\Delta f\right)=P_{\mathrm{surv}}\left(\Delta t,\Delta f\right)·𝔼\left[\sum _{m=1}^M\sum _{n=1}^N\sqrt{P_{\mathrm{qp},n_{m,}{f}_c}\left(t+\Delta t,f+\Delta f\right)}·e^{j2\pi \left(\frac{\left(f_c-f\right)d_{\mathrm{qp}},n_m\left(t\right)-\left(f_c-f-\Delta f\right)d_{\mathrm{qp}},n_m\left(t+\Delta t\right)}{c}\right)}\right]$

• • Step S704: calculating the Doppler frequency of each multipath:

$\begin{array}{c}{v}_{\mathrm{qp},n_m\left(t\right)}=-\frac{v^T}{\lambda }\cos \left({\omega }_p^T\right)-\frac{v^R}{\lambda }\cos \left({\omega }_q^R\right)\\ +\frac{{\sin }^2\left({\omega }_p^T\right){\left(v^T\right)}^{2}t}{\lambda \left(d_{n_m}^T-\cos \left({\vartheta }^T\right){\delta }_p\right)}+\frac{{\sin }^2\left({\omega }_p^R\right){\left(v^R\right)}^{2}t}{\lambda \left(d_{n_m}^R-\cos \left({\vartheta }^R\right){\delta }_q\right)}\end{array}$

• • where, δ_p/δ_q in the equation represent the antenna spacings at the transmitting end and the receiving end, ω_p^T represents an angle between the moving direction of the transmitting end and the n_mth path corresponding to the pth transmitting antenna, ω_q^R represents an angle between the moving direction of the receiving end and the n_mth path corresponding to the qth receiving antenna, and ϑ^T represents an angle between the transmitting antenna array and the n_mth path corresponding to the 1st transmitting antenna/receiving antenna; and by calculating the Doppler frequencies corresponding to all multipaths in each simulation, local Doppler spread under this simulation can be obtained, and through multiple simulations to obtain the sample mean, local Doppler spread under this type of scenario can be obtained:

$B_{qp}^{\left(2\right)}\left(t\right)=𝔼\left[\sqrt{\stackrel{\_}{\left[{v_{\mathrm{qp},n_m}\left(t\right)}^2\right]}-\stackrel{\_}{{\left[v_{\mathrm{qp},n_m}\left(t\right)\right]}^2}}\right]$

To sum up, the present invention provides a general three-dimensional geometry-based stochastic satellite channel model, which is applicable to the commonly used satellite communication bands ranging from S-band to millimeter-wave band, and, by modeling the satellite mobility based on real orbits, is suitable for the communication scenarios of LEO, MEO and GEO satellites. The model takes into account the spatial consistency of large-scale parameters and the correlation of the channel parameters: environment, frequency and satellite elevation, and supports simultaneous movement of the transmitting and receiving ends and the clusters; and it considers the birth and death of the clusters in the time domain and frequency domain, and supports the delay resolution of subpaths. In addition, the model models the ionospheric effect that is more likely to happen at the low-frequency band and the rainfall attenuation that is more likely to happen at the millimeter-wave band respectively. The present invention may be used to simulate the statistical characteristics of the channel, including effective path loss, RMS delay spread, time-frequency correlation function and local Doppler spread, and can be applied to the analysis of the effect of rain, satellite trajectory, carrier frequency and receiving-end environment on the channel, laying a foundation for the research on channel analysis, system design, etc.

It may be understood that the present invention is described with some embodiments, and those skilled in the art will appreciate that various changes or equivalent substitutions may be made to these features and embodiments without departing from the spirit and scope of the present invention. In addition, under the teaching of the present invention, these features and embodiments may be modified to adapt to a particular situation and material, without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application belong to the scope of protection of the present invention.

Claims

1. A satellite communication-oriented geometry-based stochastic channel modeling method, wherein the channel modeling method comprises the following steps: Step S1: establishing a satellite communication-oriented geometry-based stochastic channel model and a satellite channel simulation scenario corresponding to the model, and setting scenario layout parameters;Step S2: initializing the trajectory and speed of a satellite and a receiving end;Step S3: calculating spatially consistent large-scale parameters, and the effect of rainfall on the large-scale parameters;Step S4: calculating a path loss, a shadow fading, an atmospheric absorption and a rainfall attenuation;Step S5: initializing the central positions of a cluster and a scatterer, and calculating the delay, angle and power of the cluster according to the geometric position information of the transmitting and receiving ends and the scatterer, to generate a channel coefficient;Step S6: updating the large-scale and small-scale parameters according to the movement of the transmitting and receiving ends and the birth-death process of the cluster, to generate a new channel coefficient; andStep S7: deriving the statistical characteristics of the channel, and performing simulation analysis.

2. The satellite communication-oriented geometry-based stochastic channel modeling method according to claim 1, wherein the step S1 specifically comprises the following steps: Step S101: the transmitting end of the model is a satellite mobile receiving end, the receiving end is a ground mobile receiving end, the orbital altitude hsat of the satellite needs to be pre-set according to the requirements of the simulation environment, and the environment of the ground mobile receiving end also needs to be pre-set as a dense city, a city, a suburb, or a rural environment;Step S102: both the transmitting end and the receiving end of the model are a multiple-transmit multiple-receive antenna array, so the spacing and angle parameters of the antenna array need to be initialized, wherein the antenna array is a uniformly distributed linear array; the transmitting end has a total of P antennas, ApT, denotes the pth antenna, βAT denotes an azimuth angle of the transmitting-end antenna array, βET denotes an elevation angle of the transmitting-end antenna array, and a spacing between the antennas is δT; and the receiving end has a total of Q antennas, AqT denotes the qth antenna, βAR denotes an azimuth angle of the receiving-end antenna array, βER denotes an elevation angle of the receiving-end antenna array, and a spacing between the antennas is δR; andStep S103: the system parameters of the model further comprise a carrier frequency fc, and a rainfall rate set based on whether it is raining.

3. The satellite communication-oriented geometry-based stochastic channel modeling method according to claim 2, wherein the step S2 specifically comprises: Step S201: the initial position of the satellite is determined together by the orbital altitude hsat, elevation θsat and azimuth φsat; satellites are divided into GEO satellites and NGEO satellites according to whether the satellites are stationary relative to the earth; the GEO satellites are stationary relative to the ground; under the Earth's gravity, the NGEO satellites orbit the earth in an elliptical trajectory; first, an elliptical orbital plane of the satellite is determined by setting the following five parameters: a length a of the semi-major axis of the elliptical orbit; an eccentricity e of the elliptical orbit, wherein if the eccentricity e is set to 0, the elliptical orbit is reduced to a circular orbit; the angle between the orbit and the equatorial plane is an orbital plane inclination ι; the longitude of the point where the orbit crosses the equatorial plane upward is called a right ascension of the ascending node Ω; the direction of the ellipse on the orbital plane is determined by an argument of perigee ω of an angle between the orbital perigee and the ascending node; a true Anomaly υ represents an angle swept by moving along the orbit from the perigee; after determination of the orbital plane, the specific position of the satellite in the elliptical orbit can be determined via υ; the satellite trajectory is determined together by a time-varying right ascension of the ascending node Ω(t) caused by orbital perturbations and a time-varying argument of perigee ω (t) as well as a change υ(t) in the true anomaly caused by the Earth's gravity; at each moment, in the Cartesian coordinate system with the geocenter as the origin, the position coordinates (xsatysat, zsat) of the moving satellite are expressed as:

4. The satellite communication-oriented geometry-based stochastic channel modeling method according to claim 3, wherein the step S3 specifically comprises the following steps: Step S301: generating spatially consistent large-scale parameters according to the frequency, environment and satellite elevation: large-scale fading PL, shadow fading SH, delay spread DS, azimuth angular spread of arrival ASA, elevation angular spread of arrival ESA, Rician K-factor KF, and cross-polarization ratio XPR, wherein for each large-scale parameter to be calculated, the following general formula is employed:

5. The satellite communication-oriented geometry-based stochastic channel modeling method according to claim 4, wherein the step S4 specifically comprises the following steps: Step S401: calculating a free path loss PL, and modeling as a log-distance path loss model, wherein the shadow fading SF follows a log-normal distribution;Step S402: calculating the atmospheric absorption AR, setting the receiving-end height equal to the sea level height, and taking the mean annual global reference atmosphere values as the environmental parameters: temperature T, dry air atmospheric pressure p, water vapor density ρ and water vapor pressure e, wherein AG is expressed as:

6. The satellite communication-oriented geometry-based stochastic channel modeling method according to claim 5, wherein the step S5 specifically comprises: Step S501: calculating a delay τn of each cluster, a receiving-end azimuth angle φna and a receiving-end elevation angle θna using the large-scale parameters calculated in Step S301, wherein the delay of the cluster follows a unilateral exponential distribution, and an initial value thereof is calculated as:

7. The satellite communication-oriented geometry-based stochastic channel modeling method according to claim 6, wherein the step S6 specifically comprises the following steps: Step S601: introducing a survival probability Psurv(Δt, Δf) of the cluster in the time-frequency domain, with the birth and death of the cluster in the time domain and frequency domain taken into account:

8. The satellite communication-oriented geometry-based stochastic channel modeling method according to claim 7, wherein the step S7 specifically comprises: Step S701: calculating an effective path loss, which consists of three parts: free path loss PL, atmospheric absorption AG and rainfall attenuation AR, wherein the three large-scale fading components do not change rapidly with time, showing a relatively constant trend, and the sum of the three is defined as the effective path loss: PLeff=PL+AG+AR Step S702: calculating the root-mean-square (RMS) delay spread of the channel, expressed as: