Study of the Ionospheric Scintillation Radio Propagation Characteristics with Cosmic Observations

Abstract

The ionosphere has important influences on trans-ionosphere radio propagation. When signals pass through ionospheric irregularities, their amplitude and phase are often attenuated and distorted. In this work, the statistical features of scintillation observed by the Global Navigation Satellite System (GNSS) and low earth orbit (LEO) satellites are investigated with Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) data in solar cycle 24. The amplitude scintillation propagation channel is fitted by the Nakagami-m, α-μ and κ-μ models. The performance is evaluated in terms of root mean square error (RMSE), kurtosis and information entropy. The results reveal that the α-μ model achieves the best performance in all considered scintillation intensities, while the Nakagami-m model achieves better performance under severe scintillation in the GNSS-LEO propagation channels.

1. Introduction

The ionosphere refers to the upper atmosphere from 60 to 1000 km above the Earth’s surface. In the ionosphere, there are large amounts of free electrons and ions, and ionization is due to soft X-ray and far ultraviolet (UV) solar radiation as well as solar energetic particle precipitation. The ionosphere serves as a critical media for radio propagation [1]. Ionospheric irregularities occur when the electrons and ions are not uniformly distributed in the ionosphere. When radio signals pass through ionospheric irregularities, they can suffer from rapid fluctuations in amplitude and/or phase, which are named ionospheric scintillation [2]. Scintillation is impacted by solar activities and geomagnetic disturbances. It varies in accordance with geophysical locations and seasonal changes [3]. The Global Navigation Satellite System (GNSS) is influenced by scintillation. As discovered in previous studies, severe scintillation can lead to loss of signal locks, extreme ranging errors and deteriorated accuracy, integrity and availability [4,5].

There are several impact factors that are responsible for the occurrence of ionospheric irregularities [6,7]. Studies reveal that large number of ionospheric irregularities have variable patterns, such as sporadic E (Es), spread F, field-aligned irregularities (FAIs) and plasma bubbles. Ionospheric irregularities have been found in equatorial region, midlatitudes and high latitudes, due to different originations. The sporadic E layer is frequently observed from 90 to 140 km; it mostly occurs during daytime, formed in the early morning, reaching a maximum intensity at 14 LT and 20 LT for summer hemisphere [8]. The causes of E layer irregularities have been deeply investigated. Mathews attributed wind shear as important influential factor for the sporadic E features, while Farley and Simon further proposed two-stream instability (TSI) and gradient-drift instability (GDI) for E layer FAIs at equator and low latitudes. It was indicated that atmosphere gravity waves, Kelvin–Helmholtz instability, Es-layer instability and GDI were all responsible for the development of mid-latitudes irregularities in E layer [9,10,11]. Ionospheric irregularities in the F layer are more common to study, the equatorial spread F and plasma bubbles are known to be generated by Rayleigh–Taylor instability [12,13,14]. Moreover, studies reveal that atmospheric gravity waves play a critical role for the seeding of irregularities, especially during low solar activity periods. The mid-latitudes ionospheric irregularities are found usually to be associated with atmospheric gravity waves and medium-scale travelling ionospheric disturbances (MSTID) [15,16]. Its major cause is attributed to Perkins instability [7,17,18,19]. In this theory, a north–south electric field is present in addition to the east–west electric field, leading to the form of rising and falling sheets of ionization. The electro-dynamical coupling between the E layer and F layer is considered important to enhance the Perkins instability growth rate [19,20]. Geomagnetic disturbance can change the morphology of ionospheric irregularities. According to Aaron’s theory, geomagnetic storms can both excite and prohibit the formation of irregularities, depending on the storm commencement local time [21]. Previous studies on recent geomagnetic storms in solar cycle 24 indicate two dominant drivers for storm induced ionospheric irregularities, the prompt penetration electric field (PPEF) propagating from high latitudes [3,22,23] and the disturbance dynamo electric field (DDEF) formed by storm-time equatorward neutral winds and waves [6,24]. The PPEF is eastward for southward IMF Bz until after sunset (19:00 LT and later), then its direction reverses. The DDEF is westward on daytime and eastward on nighttime, opposite to the direction of PPEF.

GNSS signal propagation features under scintillation have been studied recently. Many researchers use global positioning system (GPS) signals at low latitudes to reveal the propagation channel characteristics under scintillation [25,26,27]. It was assumed that GPS signals suffered amplitude attenuations and phase jitters due to scintillation, which caused different statistical propagation distributions compared with the known Rice and Rayleigh distribution. Yeh and Liu first used the Nakagami-m model to describe the scintillation signal propagation features. The model was verified to most closely approximate the scintillation impacts on GPS signals in real experiments [1]. Wernik et al. investigated the probability distribution of ionospheric plasma intensities, revealing that the fluctuations of plasma intensities in the scintillation scale were far from the Gaussian distribution but similar to the Laplacian distribution. They also simulated the propagation of trans-ionosphere radio signals in ionospheric irregularities with two-dimensional phase screens under Gaussian and Laplacian assumptions. The relationship between non-Gaussian distributions and phase scintillation and delay time for signal propagation were discussed [28]. Moraes et al. further used real multifrequency GPS observables to study the dependency of scintillation on frequencies, with 150 nights of scintillation samples collected at the São José dos Campos site in Brazil [6]. They fitted the statistical features of GPS L1 signals under amplitude scintillation and tested the Nakagami-m, Rician and α-μ distributions. The results show the superior fitting performance of the α-μ distribution over the other two models. This work indicates that the signal features can vary under different scintillation situations, which can be represented by different combinations of α and μ parameters. Further verification was conducted with long-term ground scintillation observations in Australia, showing that the values of α and μ are related to geophysical locations. The use of entropy evaluation proved the convincing region for α-μ and Nakagami-m distributions under different scintillation intensities [29]. Based on these results, the κ-μ distribution was applied to describe the radio propagation channel under severe scintillation and compared with the α-μ distribution. It was noted that the κ-μ distribution failed to achieve a better fitting accuracy in general situations but demonstrates better approximation to reality under extreme scintillation conditions. From theoretical aspects, the α-μ, κ-μ and η-μ distributions were all proposed by Yacoub [30,31]. The α-μ distribution is considered a rewritten form of the general Gama distribution, while the κ-μ distribution is a special form of the single lateral Gaussian distribution and η-μ has potential applications for nonlinear cases. Yacoub’s contributions break the traditional constraints for signal propagation modeling under complicated situations. Further investigations under different scintillation scenarios are warranted.

Most studies on radio propagation under ionospheric scintillation have concentrated on ground-based receivers, in which the trans-ionospheric signals propagated in space-ground directions [32]. However, there is also another situation in which trans-ionospheric signals propagate from space to the receiver installed on low earth orbit (LEO) satellites, for instance, in the space-LEO directions [33,34]. GNSS signals received by Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) satellites are in such a situation. When COSMIC satellites were first launched, researchers at the University Corporation for Atmospheric Research (UCAR) tried to investigate the signal propagation paths with simulation. Launched in 2006, COSMIC satellites have been in operation for over 15 years, which has provided sufficient amounts of space-borne observations that are valuable for statistical analysis of space-LEO satellite signal propagations under ionospheric scintillation. It was discovered that the spread E irregularities led to the U shape of the received occultation signals. The U shape was closely correlated with the inclination of E irregularities [35]. Yue et al. further investigated the occultation signal features under ionospheric scintillations and statistically analyzed the loss of COSMIC satellite signal locks [36]. The COSMIC occultation signals are greatly affected by ionospheric irregularities, representing prominent scintillation characteristics. The long-term scintillation distributions on a global scale were then discussed with amplitude scintillation indices from the COSMIC received signals [37]. It was found that space-borne signal scintillation was partially correlated with the geophysical and seasonal variations in ionospheric irregularities [38].

However, there are few studies on scintillation signal distribution and models under space-LEO occultation propagation channels. To address this problem, the COSMIC scintillation observations in solar cycle 24 are processed and analyzed in this work, and several models are discussed and evaluated to fit the space-LEO satellite propagation characteristics. The COSMIC Level-2 data are used to extract the amplitude scintillation index S4 and signal-to-noise ratio (SNR) for GPS L1 frequency scintillation observations. The Nakagami-m, α-μ and κ-μ distributions are tested and fitted. The entropy function is computed to evaluate the fitting confidence levels under different scintillation intensities and in different solar activity phases. To investigate the statistical distribution of amplitude scintillation indices, the normalized SNRs of the received GPS-LEO occultation signals are fitted with the three distribution models. The fitting performances are evaluated by standard deviation, kurtosis, information entropy and fitting interval.

2. Data and Methods

2.1. The Data Presentation

The COSMIC-1 project is composed of six small satellites equipped with GPS occultation experimental devices, ionospheric photometers and three frequency beacons. It has a strong tracking ability for ionospheric scintillation. The data used in this work were collected from the GPS occultation observations provided by historical data from 2010 to 2019 in solar cycle 24, which are archived in the COSMIC Data Analysis and Archive Center (CDAAC). The COSMIC scnLv1 datasets, which contain the amplitude scintillation index (S4) and the SNR for each occultation event, are selected. Data preprocessing is implemented, including data integrity inspection, data extraction and transformation, error filtering and data classification. The amplitude scintillation index and L1 band SNR are extracted. The S4 maximum value (S4max) in 9 s along with its corresponding tangent height, latitude and longitude are recorded. The COSMIC space-borne GPS receiver does not directly measure S4; it only records the mean square root value of signal intensity fluctuation per second obtained by L1 frequency point 50 Hz carrying noise ratio data. S4 is calculated by CDAAC after the data are transmitted to the ground. Thus, errors in S4 are extracted during data preprocessing. According to [39], the amplitude of the negative fluctuation of the signal cannot exceed the signal strength, and the amplitude of the positive fluctuation cannot exceed the mean value of the signal for a long time. For the convincing aspect of the analysis, the range of S4 is set between 0.3 and 1.5, eliminating S4 larger than 1.5. The threshold of occultation scintillation is set to 0.3 since a scintillation intensity below 0.3 is of little significance to this work. The S4 values are further normalized to make them suitable for statistical modeling and intensity comparisons. After normalization, the scintillation intensity is then classified into different levels, i.e., very weak, weak, moderate and strong (

Table 1. Ionospheric scintillation intensity levels.

S4	Ionospheric Scintillation Intensity
<0.3	Very weak
0.3–0.5	Weak
0.5–0.7	Moderate
0.7–1	Strong

).

In previous studies, the signal amplitude was derived by the carrier-to-noise ratio (CNR). Because it is difficult to determine the CNR and bandwidths for COSMIC data, the SNR is used and normalized to extract the distribution of the signal amplitude. To exclude noise interference, an SNR higher than a threshold of 15 dB is considered and further processed. When the amplitude distribution in relation to S4 is computed, the SNR is classified by the given S4 with ±0.005 intervals to extract sufficient SNR samples. The narrower interval of S4 variation is more conducive to improving the accuracy of amplitude statistical fittings.

Ionospheric scintillation depends on the geophysical locations. The most representative scintillation events are concentrated around geomagnetic latitudes of ±20°, where the equatorial ionization anomaly (EIA) exists. The long-distance lateral transmission of occultation signals makes the scintillation observations subject to the integration of electrons with uneven paths, leading to difficulties in judging the precise location of the irregularities that are responsible for scintillation. In this work, the observed S4 from COSMIC data are grouped by tangent geographical latitudes into 0–30°, 30–60° and 60–90°, representing low latitudes, including the EIA region, middle latitudes and high latitudes, respectively. To characterize the scintillation in vertical aspects, the tangent height is divided into two groups, 100–150 km and 150–500 km, to distinguish the scintillation characteristics between the E and F layers, respectively. The COSMIC observation satellites at different latitudes and different heights are distributed heterogeneously. Thus, simply using the amount of data to analyze the occurrence times of scintillations of different intensities is not appropriate. The occurrence probability of different intensities of scintillation is calculated and used to analyze the latitude and height distribution of scintillation characteristics.

2.2. Ionospheric Scintillation Intensity

The amplitude scintillation index is defined as [40]

$$S4=\sqrt{\frac{⟨I^2⟩-{⟨I⟩}^2}{{⟨I⟩}^2}}$$

(1)

where S4 indicates the amplitude scintillation intensity. $I=R^2$ represents the strength of the received signal. $R$ is the signal amplitude. $⟨·⟩$ indicates the time average. Amplitude scintillation occurs at a higher probability at low latitudes, while it seldom occurs at high latitudes.

The phase scintillation index is given as

$${\sigma }_{\varphi }=\sqrt{{\varphi }^2-{\varphi }^2}$$

(2)

where $\varphi$ is the tracked signal phase after detrending. Severe phase scintillation mainly occurs at high latitudes and polar regions.

There are different distributions to describe the scintillation signal propagations, while the most traditional distribution is the Nakagami-m distribution, given as [41]

$$f_R\left(r\right)=\frac{2m^m{r}^{2m-1}}{\Gamma \left(m\right){\Omega }^m}{e}^{-m{r}^2/\Omega }$$

(3)

where r is the unit vector of the received signal amplitude vector $R$. $\Gamma \left(·\right)$ indicates the Gama function. $m=1/{\left(S_4\right)}^2$ represents the attenuation. $\Omega =E\left\{r^2\right\}$ is the average power of received signals. After the normalization, Equation (3) is rewritten as

$$f_R\left(r\right)=\frac{2m^m{r}^{2m-1}}{\Gamma \left(m\right)}{e}^{-m{r}^2}$$

(4)

The α-μ distribution is considered to be a rewritten form of the generalized Gama distribution. It is used to describe the distribution of the GPS scintillation signals at low latitudes in America, where amplitude scintillation is the most prominent [25]. This model assumes the signal as a composed of clusters of multipath waves propagating in a nonhomogeneous environment. This distribution has two degrees of freedom and is thus more suitable for describing the amplitude distribution under scintillation. The α-μ distribution is given by [42]

$$f\left(r\right)=\frac{\alpha r^{\alpha m-1}}{{\xi }^{\alpha m/2}\Gamma \left(\mu \right)}exp\left(-\frac{r^{\alpha }}{{\xi }^{\alpha /2}}\right)$$

(5)

where $\xi$ is calculated by $\xi =\frac{\Gamma \left(\mu \right)}{\Gamma \left(\mu +2/\alpha \right)}$. The α, μ parameters are estimated by [42]

$$\frac{E^2\left(R^{\beta }\right)}{E\left(R^{2\beta }\right)-E^2\left(R^{\beta }\right)}=\frac{\Gamma \left(\mu \right)\Gamma \left(\mu +4/\alpha \right)-{\Gamma }^2\left(\mu +2/\alpha \right)}{{\Gamma }^2\left(\mu +2/\alpha \right)}$$

(6)

where $\beta$ is a parameter to be determined. As there are two unknowns, assuming $\beta =3$ and $\beta =4$, Equation (6) can be rewritten as two equations, and $\alpha$ and $\mu$ can be calculated.

It should be noted that the same $S_4$ can correspond to different α parameters. A larger value of α indicates that the signal is more attenuated. When $\alpha =2$ and $\mu =m$, the α-μ distribution equals the Nakagami-m distribution. Thus, Nakagami-m is a special case of α-μ distributions.

The α-μ distribution is a kind of envelope distribution that uses two parameters, α and μ, to establish a generalized fading model considering multipath effects. It assumes that the signals propagating in a nonhomogeneous environment have similar delay times, with delay-time spreads of different clusters being relatively large. The scattered signals diverge over certain distances, but all of the scattered signals have the same impact on the phase, with the same power. Then, the envelope can be abstracted as a nonlinear function of the sum of multipath scatters. The parameter α indicates the extent of nonlinearity, while the parameter μ is related to the number of multipath scatters. The α-μ distribution takes both the nonlinearity effect and wave cluster effect into account. With increasing amplitude scintillation index $S_4$, the peak value of the α-μ distribution decreased, and with a corresponding smaller amplitude or SNR, the curve diverged. Under the same $S_4$, a larger α will decrease the peak value of the α-μ distribution. The width and thickness of the curve tail also increased. The parameters α and μ have hyperbolic-like features. The α increases with the decrease of μ.

The κ-μ distribution model is developed based on the α-μ distribution hypothesis that the signal is a combination of multipath signals passing through heterogeneous propagation media [27,31]. The model is defined as [31]

$$f_R\left(r\right)=\frac{2\mu {\left(1+\kappa \right)}^{\left(\frac{\mu +1}{2}\right)}}{{\kappa }^{\left(\frac{\mu -1}{2}\right)}{e}^{\kappa \mu }}{r}^{\mu }{e}^{\left[-\mu \left(1+\kappa \right)r^2\right]I_{\mu -1}\left[2\mu \sqrt{\kappa \left(1+\kappa \right)}r\right]}$$

(7)

where $I_{\mu -1}$ represents the modified first kind Bessel function of order $\mu -1$. $\kappa$ represents the ratio between the dominant components versus the power of the scattered wave, which can be calculated by

$$\kappa =\frac{\sqrt{2}\left[E\left(R^4\right)-1\right]}{\sqrt{2E^2\left(R^4\right)-E\left(R^4\right)-E\left(R^6\right)}}-2$$

(8)

Similar to the α-μ distribution, the parameters κ and μ are related to the S4 index via

$$S_4=\sqrt{\frac{1+2\kappa }{\mu \left(1+{\kappa }^2\right)}}$$

(9)

2.3. Parzen Window

The Parzen window is referred to as the kernel density estimation method to realize the nonparameter test. Suppose that the samples follow an unknown probability density function while lying in D dimension space $R$. The probability of each sample falling into the $R$ space is defined as

$$P={{\displaystyle \int }}_{}^{}p\left(x\right)dx$$

(10)

where $p\left(x\right)$ is the probability density function. Then, we employ k to represent the number of samples in space $R$. Equation (8) is further written as $p\left(x\right)\approx \frac{k}{nV}$, where n is the number of samples and $V$ is the volume in space $R$. The window function is defined as

$$H\left(\mu \right)=\left\{\begin{array}{c}1,\left|\mu \right|\le 0.5\\ 0,otherwise\end{array}\right.$$

(11)

$k$ is then calculated as

$$k_n={\displaystyle \sum }_{i=1}^{n}H\left(\frac{x-x_i}{h_n}\right)$$

(12)

where $h_n$ is the length of space $R$. The probability density is then solved as

$$P_n\left(x\right)=\frac{1}{nh}{\displaystyle \sum }_{i=1}^{n}H\left(\frac{x-x_i}{h_n}\right)$$

(13)

This method is often used to estimate the probability density function of samples. However, the estimated density equation is not smooth enough, and the estimation result is affected by the group distance. In this work, the probability density function and phase distribution of scintillation under different intensities are estimated by the Parzen window method.

2.4. Root Mean Square Error (RMSE)

In the theory of mathematical statistics, residuals and variances are often used to evaluate the degree of fitting and regression of the probability density function (PDF). In this study, RMSE is analyzed to systematically evaluate the difference between the measured and theoretical distributions.

The RMSE is given as

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(14)

where $y_i$ is the observation, ${\widehat{y}}_i$ is the estimation, and RMSE is a widely used metric to evaluate the difference between the estimation and the real value.

2.5. Kurtosis and Information Entropy

Kurtosis is introduced to explain the contour and non-Gaussian properties of the probability density distribution. It takes the normal distribution as the standard to determine the tail distribution of a specific PDF. A greater value of kurtosis indicates a thicker tail and a wider range of the curve. Kurtosis is used to study the properties of the signal amplitude distribution under different scintillation intensities. For a fourth-order random variable, the kurtosis is given as [43]

$$C_k=\frac{E{\left(x-E\left(x\right)\right)}^4}{{\left[E{\left(x-E\left(x\right)\right)}^2\right]}^2}-3$$

(15)

where $x$ is sample values.

Information entropy was originally used in the field of informatics to measure the disorder of information and reflect the amount of information at the same time. Currently, it is used to quantitatively describe the degree of chaos and uncertainty of the system. The higher the value is, the higher the disorder degree of the system. Information entropy is expressed as

$$H=-{\displaystyle \sum }_{i=1}^m{P}_i{\log }_2{P}_i$$

(16)

where $P_i$ indicates the probability. The equation fails when $P=0$. Since ionospheric scintillation is a random event with high disorder and uncertainty, information entropy is used here to discuss the confidence level of the amplitude distribution under different scintillation intensities. A similar usage of information entropy in scintillation studies is referred to in Guo et al. 2017 [29].

2.6. Fitting Interval Analysis

The interval of samples is also considered in the evaluation of fitting performance. Due to the errors in the data or induced during the data processing, the “mean ± standard deviation” is commonly used to evaluate the fitting performance. In this study, the value of “mean ± standard deviation” was adopted to count the interval length of the standard curve falling into this interval.

The proposed framework of this analysis is described in Figure 1. The COSMIC occultation data are first downloaded, and then the right scintillation index is extracted. The data are classified to different levels. The Parson window method is then adopted to calculate the probability density function. The confidence interval is also determined. Finally, the characteristics and fitting performance of the measured curve and the standard probability density distribution curve are compared and discussed by using the statistical metrics introduced above.

3. Experimental Results

3.1. Temporal Statistical Feature of Amplitude Scintillation Index S4

Figure 2 shows the histogram of yearly amplitude scintillation occurrence of different intensities in layers E and F. The probability of scintillation for each intensity in layer E changes slightly. Overall, weak scintillation accounts for the largest proportion (approximately 69%), while moderate scintillation ranks second and accounts for approximately 20%. The strong scintillation intensity is only approximately 11%. In the F layer, the weak scintillation intensity occurrence generally increases year by year, accounting for approximately 70% on average. The moderate and strong scintillation occurrences decrease with decreasing solar cycle, taking up 20% and 10%, respectively. In both the E and F layers, the scintillation occurrences decrease with the increasing in intensity. This is consistent with the conclusions of previous studies [25,26,27] that the probability of scintillation decreases gradually with increasing scintillation intensity.

Based on Figure 2, the scintillation occurrence feature observed in the ionospheric E and F layers was summarized; it was found that the occurrence of weak scintillation in the E layer is lower than that in the F layer, indicating that weak scintillation occurs relatively less frequently in E layer. By contrast, for moderate and strong scintillations, the occurrence in the E layer is greater than that in the F layer. In general, with the increase in scintillation intensity, the amplitude scintillation occurrence decreased gradually.

Figure 3 shows the occurrences of scintillation in the E and F layers at different latitudes from 2010 to 2019. The occurrence of different intensities of scintillation at low, middle and high latitudes in the E layer is almost the same. In the F layer, the weak scintillation occurrences at different latitudes show a similar trend, with a concave trend in 2012, which is the solar maximum in this solar cycle. For moderate and strong scintillation, similar trends are noticed. They first increased from 2010 to 2014 and then gradually decreased corresponding to the variation in solar activities.

3.2. Statistical Analysis for SNR

Figure 4 shows the probability density function of the normalized SNR when the normalized S4 values are 0.5, 0.7 and 0.9. In statistical analysis, the Parzen window method is used to fit the measured curve. Here, the fitting result is called the evaluated curve, which is applied to calculate the Nakagami-m, α-μ and κ-μ distribution models. All the data from solar maximum year 2014 and solar minimum year 2019 are selected for analyses. When S4 is 0.5, with moderate scintillation, the fitting curves of the Nakagami-m distribution and κ-μ distribution are nearly the same. On the other hand, the α-μ distribution better approximates the measured SNR in 2019, when the solar activity is in the low phase. When S4 is 0.7, the fitting curves of the Nakagami-m distribution and κ-μ distribution begin to diverge, and the α-μ distribution still outperforms the other two distributions in approximating the measured SNR in both 2014 and 2019. When S4 is 0.9, the fitting performances of these three distributions are the most similar in 2014, with the α-μ distribution achieving the best fitting performance. For 2019, the α-μ distribution loses efficiency, and the κ-μ distribution outperforms the Nakagami-m distribution.

Figure 4 shows that the distribution of the normalized SNR demonstrated non-Gaussian distribution features [35]. The protruding part shows another peak of the distribution curve. With the increase in normalized S4, the protrusion decreases gradually. When the normalized S4 reaches approximately 0.7, the protrusion becomes mild.

The ten-year COSMIC occultation data are studied, and the corresponding SNR distributions are shown in Figure 5. The Nakagami-m and κ-μ distribution curves are quite close to each other when S4 ranges from 0.4 to 0.6, showing “high and thin” shapes. The α-μ distribution fits the measured SNR better when $S4\ge 0.6$. With increasing scintillation intensity, the fitting performance of the Nakagami-m and κ-μ distributions starts to diverge. When $S4$ is 0.9, the fitting performances of the Nakagami-m and α-μ distributions are very similar, and when $S4$ is 1.0, the Nakagami-m distribution outperforms the other two distributions. It is preliminarily considered that the fitting performance of the α-μ distribution curve is the best when $S4\ge 0.6$, which is consistent with the results from ground-based observations [29].

Table 2. Calculated values of parameters α and μ under different scintillation intensities.

S4	α	μ
0.3	5.3467	0.239
0.4	4.2748	0.3481
0.5	3.788	0.3846
0.6	3.1258	0.514
0.7	2.5877	0.8942
0.8	2.3442	1.0059
0.9	1.9489	1.3512
1.0	1.4625	2.1049

shows the calculated values of the parameters α and μ under different scintillation intensities assuming β = 3. As shown by Moraes et al., the relationship of α and S4 can be estimated [41]. By analyzing the α parameters in the table, it is noticed that the α parameter decreases with the increase in S4. This phenomenon also occurs in the variation of the α parameter calculated in ground-based observations. When the normalized amplitude scintillation index is equal to 0.9, α is 1.9489, close to 2. When α = 2, the Nakagami-m distribution function and α-μ distribution function are equal. Therefore, the fitting performance of these two standard probability density distribution curves should be close [43].

3.3. Goodness of Fit for Different Distributions

In this part, residuals of the estimation deviation are analyzed. Figure 6 shows the standard deviation of the normalized amplitude scintillation index for three distributions under different intensities. This result indicates that a better fitting affect is achieved the smaller the values of the standard deviation. The difference between the Nakagami-m distribution curve and measured probability density curve decreases with increasing normalized S4, indicating that the fitting performance of the Nakagami-m distribution curve improves with increasing normalized S4. The results are different from those obtained from ground-based observations. From the ground-based results, when normalized $S4\le 0.7$, the standard deviation between the Nakagami-m distribution function curve and the measured probability density curve decreases with increasing S4, showing that the fitting performance of the Nakagami-m distribution function is better with increasing S4 under weak and medium scintillation [29]. When the normalized S4 is greater than 0.7, the root mean square error of the Nakagami-m distribution has an upward trend in the ground-based results. This reflects the difference between the analysis results of radio occultation data and ground data. The standard deviation of the α-μ distribution is generally low, with more slight fluctuation, indicating that the fitting effect of the α-μ distribution is relatively better as a whole.

When the normalized S4 is equal to 0.9, the standard deviation of the two analytical distribution curves basically overlaps. This is because when the normalized S4 is equal to 0.9, the value of α approaches 2, and the two probability density distribution curves are basically equal, so the standard deviation is also similar. When the normalized S4 is greater than 0.9, the standard deviation of the Nakagami-m distribution is less than that of the α-μ distribution, indicating that the fitting performance of the Nakagami-m distribution exceeds that of the α-μ distribution when the ionospheric scintillation intensity is extremely severe.

The results of the κ-μ distribution demonstrate a similar variation trend of the standard deviation as that of the Nakagami-m distribution, and intersection occurs when the normalized S4 is 0.5. It also shows that the fitting performance of the κ-μ distribution improves with increasing normalized S4, but the fitting performance is not as good as that of the Nakagami-m distribution when the scintillation intensity becomes extremely severe.

The kurtosis and information entropy analysis methods are used to analyze the fitting performance and evaluate the actual fitted probability density curve.

Kurtosis reflects the nonuniformity of the curve tail distribution. The greater the kurtosis value, the wider the tail. When studying the statistical distribution of SNR under a certain S4, it is observed that the probability density distribution does not fit the normal distribution. In some actual curves, the two ends of the normalized amplitude do not decrease monotonically, but there are abnormal peaks in some places. The kurtosis difference (absolute value) between three standard distribution functions and the measured curve under different normalized amplitude scintillation indices is calculated as shown in Figure 7.

For the Nakagami-m distribution, with the increase of the normalized S4, the difference of kurtosis decreases, indicating that the fitting performance of Nakagami-m distribution in this range is getting better. Meanwhile, for the α-μ distribution, the kurtosis difference is generally small regardless of the change in the normalized S4, and the value is mostly below 1. In the range of normalized S4 from 0.8 to 0.9, the kurtosis difference between the Nakagami-m and α-μ models intersects, and then the kurtosis of Nakagami-m becomes less than that of α-μ, indicating that when $S4>0.9$, the Nakagami-m distribution function is better than the standard α-μ distribution. For the kurtosis difference curve of the κ-μ distribution, when the normalized S4 is less than 0.5, it is below the kurtosis of the Nakagami-m distribution. When the normalized S4 is greater than 0.5, it exceeds the kurtosis of the Nakagami-m distribution.

Information entropy reflects disorder and uncertainty of specific results. In this study, information entropy is used to evaluate the disorder of the distribution curve. Figure 8 shows the information entropy difference (absolute value) between the measured probability density distribution and the theoretical distribution under different normalized S4. With the increase in normalized S4, the information entropy difference of the Nakagami-m distribution generally exhibits a downward trend, and there are still some fluctuations. In the case of weak ionospheric scintillation, the disorder between the Nakagami-m distribution and the actual signal-to-noise ratio distribution is quite different. When S4 increases to 1.0, the difference in information entropy reaches the minimum, indicating that the signal-to-noise ratio distribution is affected by intensity modulation, and the disorder difference has reached a low level. In the information entropy difference curve of the α-μ distribution, it is found that the size of the information entropy difference is generally small. The two curves also overlap when the normalized S4 is 0.9, indicating that the results of the goodness of fit evaluated by information entropy matches the results in Figure 5. For the κ-μ distribution, the information entropy difference overlaps with the Nakagami-m distribution when S4 is 0.5, which is consistent with the results in Figure 5.

To study the performance of the three standard distributions in different years in solar cycle 24,

Table 3. RMSE, information entropy and kurtosis for the three distributions in 2014 and 2019 ($S4=0.7$).

Year	RMSE			Information Entropy			Kurtosis
	Nakagami-m	α-μ	κ-μ	Nakagami-m	α-μ	κ-μ	Nakagami-m	α-μ	κ-μ
2014	0.1405	0.0708	0.2043	0.5060	0.1365	0.727	0.7146	0.0245	0.0381
2019	0.1459	0.0993	0.2021	0.6609	0.0797	0.8777	0.1145	0.6141	0.2618

shows the standard deviation, information entropy and kurtosis difference between the three distribution curves and the measured probability density distribution curves in 2014 and 2019 when the normalized S4 is set to 0.7.

When normalized S4 is 0.7, the standard deviation and information entropy difference between the Nakagami-m distribution and the measured normalized SNR probability density curve in 2014 (solar maximum year) are slightly less than the corresponding values in 2019 (solar minimum year), while the kurtosis difference is larger than the corresponding values in 2019. This demonstrates that, from the perspective of standard deviation and information entropy, the α-μ distribution fits better in high solar activity years. When normalized S4 is 0.7, the standard deviation and kurtosis difference between the α-μ distribution in 2014 and the measured normalized signal-to-noise ratio probability density curve are less than the corresponding values in 2019, while the information entropy difference is larger than the corresponding values in 2019. When the normalized S4 is 0.7 in 2014, the distribution for κ-μ is close to the standard deviation of the measured normalized SNR probability density curve. The information entropy and kurtosis differences are less than the corresponding values.

When the normalized S4 is 0.7, the standard deviation, information entropy difference and kurtosis difference between the α-μ distribution and the measured normalized SNR probability density curve are less than the corresponding values of the other two distributions in 2014. The fitting effect of the κ-μ distribution is the least ideal, but the corresponding kurtosis value is close to that of the α-μ distribution, achieving a better fitting of the tail. When the normalized S4 is 0.7, the standard deviation and information entropy difference between the α-μ distribution and the measured normalized SNR probability density curve in 2019 are less than the corresponding values of other distributions, while the kurtosis difference is the largest of all. From the perspective of kurtosis, the Nakagami-m distribution has a better fitting performance in 2019, while from the perspective of standard deviation and information entropy, the α-μ distribution has a better fitting performance in 2019. In both 2014 and 2019, the standard deviation and information entropy of the κ-μ distribution are the largest. This means that from the perspective of standard deviation and information entropy, the κ-μ distribution fits the worst of the three. The kurtosis indicates that the κ-μ distribution achieves a better fitting of the tail.

A further investigation was conducted to divide the results into two groups: an ascending solar phase from 2012 to 2015, and a descending solar phase from 2016 to 2019. At this time, the normalized S4 ranged from 0.3 to 1.0. Figure 9 shows the difference between the three models and the measured probability density curve. Slight differences in RMSE, information entropy and kurtosis are present, which means that for a general situation with large samples, the influence of solar activity on a single model turns to be very slight. The α-μ model outperforms in all considered aspects, then followed by the Nakagami-m and κ-μ distributions.

The fitting confidence interval is evaluated, and Figure 10 shows the result of $S4=0.7$. Both α-μ and Nakagami-m are between the upper and lower limits, while for κ-μ, a part of the probability density function curve exceeds the fitting confidence boundary. This indicates that the κ-μ model cannot achieve good performance when the amplitude scintillation is strong (but not very severe), consistent with the above results.

By comprehensively analyzing the standard deviation, kurtosis and information entropy, the fitting performance of the measured probability density distribution curve can be best fitted with the α-μ distribution. The “mean ± standard deviation” is used to estimate one standard deviation for the estimation curve. The interval length and probability of the three distributions falling into the range are counted to find the best fitting interval.

Table 4. Intervals falling into the “mean ± standard deviation” and the corresponding interval probability of the fitting performance under different scintillation intensities.

S4	Intervals in Nakagami-m Distribution	Interval Probability in Nakagami-m Distribution	Intervals in α-μ Distribution	Interval Probability in α-μ Distribution	Intervals in κ-μ Distribution	Interval Probability in κ-μ Distribution
0.3	(0,0.32) (0.73,0.8) (1.23,1.35) (1.82,3.63)	63.91%	(0,3.63)	100%	(0,0.32) (0.67,0.77) (1.25,1.43) (1.82,3.63)	66.39%
0.4	(0–0.29) (0.65,0.76) (1.24,1.41) (1.78,4.07)	70.27%	(0,4.07)	100%	(0,0.29) (0.61,0.75) (1.25,1.45) (1.78,4.07)	71.74%
0.5	(0,0.23) (0.56,0.72) (1.26,1.49) (1.81,4.22)	71.8%	(0,4.22)	100%	(0,0.23) (0.55,0.73) (1.26,1.5) (1.81,4.22)	72.51%
0.6	(0.4,0.64) (1.25,1.52) (1.86,4.27)	68.38%	(0,4.27)	100%	(0.44,0.67) (1.26,1.5) (1.88,4.27)	66.98%
0.7	(0,0.8) (1.1,3.65)	91.78%	(0,3.65)	100%	(0,0.76) (1.17,3.65)	88.77%
0.8	(0,3.54)	100%	(0,3.54)	100%	(0,0.76) (1.16,3.54)	88.7%
0.9	(0,3.76)	100%	(0,3.76)	100%	(0,0.77) (1.16,3.76)	89.63%
1.0	(0,3.72)	100%	(0,3.72)	100%	(0,0.75) (1.13,3.72)	89.78%

lists the intervals falling into the “mean ± standard deviation” of the fitting performance of the three standard distributions and the corresponding interval probability under different scintillation indices. The data are comprehensively processed from 2010 to 2019.

The fitting range of the Nakagami-m distribution becomes larger with increasing the normalized S4 and approaches 100% when $S4\ge 0.8$. The fitting range of the α-μ distribution satisfies 100% goodness of fit in the whole range, while the fitting range of the κ-μ distribution has a similar trend to the Nakagami-m distribution but cannot meet 100% goodness of fit in the whole range.

Figure 11 shows the profile of the skewness and kurtosis. In other studies, it has been found that the profile is parabolic centered on zero, indicating that skewness and kurtosis have an empirical quadratic relationship [44]. In relevant ionospheric studies, this is considered an indication of the local dynamics of the plasma within the scattering ionospheric layer [45]. In this work, we also investigate the skewness and kurtosis profile using the data from 2010 to 2019 of all scintillation intensities, and the profile tends to collapse to the parabolic line, which is coincident with previous studies.

4. Discussion

From the above experiments, it is considered that the overall fitting performance of the α-μ distribution is the best among the three. The α-μ distribution adopts two parameters in fitting the signal SNR, increasing the degree of freedom. The standard Nakagami-m distribution only uses one parameter, but it achieves good fitting performance when $S4\ge 0.8$. Although the κ-μ distribution also uses two parameters, its performance is limited, and it only achieves better performance when $S4\le 0.5$. The α-μ distribution fits well in the whole range of the normalized S4 with a 100% confidence level. The Nakagami-m distribution continues to improve with increasing normalized S4, indicating that the Nakagami-m distribution is more suitable under strong scintillation intensities. The fitting performance was also evaluated in high and low solar activity years, which shows that the results were influenced by solar activity, but not by much. A probable reason is that more moderate and strong scintillation occurred in the high solar activity year (2014 for example), as indicated by Figure 3.

When using the standard deviation, information entropy and kurtosis for fitting evaluation, there are consistent conclusions. The α-μ distribution has the best performance for all scintillation intensity ranges, the Nakagami-m distribution gradually achieves better performance with increasing scintillation intensity, and the κ-μ distribution outperforms the Nakagami-m distribution under weak scintillations.

The amplitude scintillation distribution for the GPS-LEO propagation channel shows different features from that for the GPS-ground propagation channel, as revealed by Figure 5 and Figure 6. For the Nakagami-m distribution, the fitting performance improves under strong and severe scintillation for the GPS-LEO propagation channel, while for the GPS-ground propagation channel, the Nakagami-m distribution achieves good performance for moderate scintillation. The possible reason is that for the GPS-LEO propagation channel, the receiving signal travels a bending path, with more interaction with the ionosphere than the GPS-ground propagation condition.

5. Conclusions

In this work, the global ionospheric scintillation feature for radio propagation is investigated with the COSMIC occultation data from 2010 to 2019. The statistical characteristics of the amplitude scintillation index $S4$ and SNR are analyzed, and the radio propagation channel feature is extracted. It was found that the probability for weak scintillation accounts for the largest proportion, with approximately 69%, while the moderate and strong scintillation accounts for approximately 20%, and 11% respectively.

The SNR distribution curves under different normalized $S4$ are analyzed in each year. When the normalized S4 is small, there is a bimodal phenomenon in the measured distribution curve, which is probably caused by different factors affecting ionospheric scintillation, such as atmospheric multipath refraction. The Nakagami-m, α-μ and κ-μ distributions are used to fit the amplitude scintillation feature. From the results, it is considered that the overall fitting performance of the α-μ distribution curve is the best. When the normalized $S4\ge 0.8$, the Nakagami-m distribution function achieves better performance. The κ-μ distribution is suitable for weak scintillation cases. The α and μ values under different normalized S4 are calculated and analyzed by confidence intervals.

To evaluate the fitting performance of the three distribution curves, the standard deviation, kurtosis, information entropy and confidence interval are used. Through the analysis of standard deviation, the following are found: (a) With the increase of the normalized S4, the standard deviation of difference between the Nakagami-m distribution and the measured curve decreases, indicating that the fitting performance of the Nakagami-m distribution becomes better with increasing normalized S4. (b) The distribution curve fitted by the α-μ distribution has a generally smaller standard deviation difference than the other two distributions. (c) The standard deviation of the difference for the κ-μ distribution shares the similar variation trend as that of the Nakagami-m distribution, while the fitting performance is not as good as that of the Nakagami-m distribution. In the Nakagami-m distribution, the kurtosis difference decreases with increasing normalized S4, indicating that the fitting performance of the Nakagami-m distribution in this range improves. In the α-μ distribution, the kurtosis difference is generally small regardless of the change in normalized S4. The κ-μ distribution is outstanding among the three distributions to fit the low scintillation intensity case. A similar trend was noticed for the information entropy evaluation.

The fitting performance of the three distribution models is analyzed in different solar activity years. The fitting interval analysis was discussed to evaluate the fitting performance. The fitting interval length of the Nakagami-m distribution increases with increasing normalized S4; α-μ distribution meets the 100% confidence level in the whole range; κ-μ distribution has a similar trend as that of the Nakagami-m distribution, while it does not reach 100% in the whole range. The skewness–kurtosis profile extracted in this work shows similar features as revealed by a previous study, indicating that the amplitude scintillation propagation channel is probably related to the nonlinearity of ionospheric irregularity and potential influential factors of the ionosphere.