# Estimation of Scintillation Indices: A Novel Approach Based on Local Kernel Regression Methods

Authors: Mohammed Ouassou; Oddgeir Kristiansen; Jon G. O. Gjevestad; Knut Stanley Jacobsen; Yngvild L. Andalsvik
Journal: International Journal of Navigation and Observation
 (2016)

---

## Abstract

We present a comparative study of computational methods for estimation of ionospheric scintillation indices. First, we review the conventional approaches based on Fourier transformation and low-pass/high-pass frequency filtration. Next, we introduce a novel method based on nonparametric local regression with bias Corrected Akaike Information Criteria (AICC). All methods are then applied to data from the Norwegian Regional Ionospheric Scintillation Network (NRISN), which is shown to be dominated by phase scintillation and not amplitude scintillation. We find that all methods provide highly correlated results, demonstrating the validity of the new approach to this problem. All methods are shown to be very sensitive to filter characteristics and the averaging interval. Finally, we find that the new method is more robust to discontinuous phase observations than conventional methods.


---

## Body

## 1. Introduction

Global Navigation Satellite Systems (GNSS) are based on satellite signals being transmitted to receivers on the ground. These signals must pass through the ionosphere on the way, where the variation in electron density causes undesirable fluctuations in the observed signal. This distortion, known asscintillation, can affect both the amplitude and phase of GNSS signals. The relevant electron density variations range from decameters to kilometers in size [1–3].In the auroral zones, it is much more common for phase scintillation than amplitude scintillation to be the dominant effect (e.g., [4–6]). Scintillation is a large challenge for navigational systems, as it can disturb not only single-receiver systems but also networked systems [3, 4, 7–9]. In a worst-case scenario, the user receiver can end up with a loss of signal lock, which causes discontinuities in the phase measurements. This is known as acycle-slip.A good understanding of the ionosphere morphology can aid the development of suitable mitigation algorithms, that is, algorithms that assign weights to observations based on the scintillation distortion of each satellite link.To be more specific, a more realistic stochastic model for GNSS observables would have to take into account the variance caused by scintillation. This would avoid biased solutions. Today, the stochastic model includes the following: correlation among observations [10]; satellite elevation dependency modeled by exponential function [11]; temporal and cross-correlations [12]; and multipath detection and monitoring [13]. That is, a suitable robust weighting algorithm that reduces the influence of the satellite exposed by scintillation will look similar to the one proposed by Eueler and will enhance the ability to resolve the carrier-phase ambiguity and to improve the stochastic model for GNSS processes.This is most commonly gauged by thephase scintillation index
σ
ϕ, which is simply defined as the standard deviation of the detrended carrier phase ϕ over some period of time. A critical step in the calculation of σ
ϕ is a frequency filter, which is used to separate the high-frequency ionospheric distortions from the low-frequency distortions due to, for example, multipath interference. The scintillation indices are usually calculated for one-minute intervals and then processed with a Butterworth sixth-order high-pass filter [14–23]. However, the estimated value of σ
ϕ can be highly sensitive to the cutoff frequency of the filter. While the most common cutoff frequency is 0.1 Hz, it has been shown that this is suboptimal at high latitudes and that higher values like 0.3 Hz might yield better results [24]. The variation in the carrier phase of GNSS signals is also commonly quantified using the rate of total electron content index (ROTI) [25–27]. Several other indices for quantifying phase variations have also been proposed [24, 28–32]. It is worth noting that all of these indices are directly related to the scintillation signal variance.Some scintillation indices use methods based on wavelet transforms for the filtration. A major advantage of wavelet filters is that they manage to preserve local features in the signal, thus avoiding misinterpretation that can occur when using standard filter approaches [24, 33–36]. While such wavelet filters may yield better results than conventional filters, the computational load can become very high, especially when processing data with high sampling rates (50+ Hz). The improvement that can be achieved by substituting a wavelet filter for the conventional filters must therefore be weighed against the computational load required, especially if the algorithm is intended for real-time applications.The primary focus of this paper is on nonparametric local regression with bias Corrected Akaike Information Criteria as a viable alternative for the computation of reliable scintillation indices. In addition, we perform a statistical analysis of the scintillation indices and show that they can be well described by the Nakagami and Frechet distributions. The paper is organized as follows.Section2: it includes detailed description of the data sets used in this paper; Section 3: it deals with digital filter construction of low-pass filter algorithms for amplitude scintillation index computation and high-pass filter algorithms for phase scintillation computation. Filter types discussed are the Chebyshev, Butterworth, and Elliptic filters. These algorithms are used to validate the new approach; Section 4: it tackles nonparametric regression with kernel smoothing; Section 5: this section presents how the scintillation indices are computed by the conventional methods and the new approach; Section 6: it shows statistical analysis of scintillation indices, including distributions, sensitivity, and correlation analysis of implemented algorithms; Section 7: it presents generalization of the preceding sections; Section 8: it includes conclusion and applications of the results.


## 2. Test Data

The data used for the investigation herein was obtained from three different NRISN observation sites in the northernmost parts of Norway. These sites are highlighted with an italic typeface in Table1 and encircled in Figure 1. The observational data is taken from day 50 of year 2014, which according to the I95 index corresponds to relatively high ionospheric activity. (TheI95 index  [37] is used for the classification of ionospheric activity by CPOS software.) For more information about the ionospheric activity levels, see Figures 15 and 16. The NRISN sites are equipped by dual-frequency Septentrio PolaRxS receivers tracking GPS and GLONASS L1 and L2 frequencies with sampling rate of 100 Hz and are capable of tracking up to 12 satellites for each system. The baselines vary between 176 and 1988 km, and the height difference between sites is around 400m. To give a full picture of the baseline lengths and the height difference between sites above the ellipsoid WGS84, Tables 2 and 3 are provided.Table 1
List of NRISN observation sites.
Site
4-char
ID
Receiver type
Antenna type
Tromso
TRO2
1001
Septentrio PolaRxS
MCGNSS
Septentrio
Vega
VEG2
1002
Septentrio PolaRxS
PolaNt_MC
Ny-Ålesund
NYA2
1003
Septentrio PolaRxS
JAV_RING-ANT_G3T
Hofn
HOF2
1010
Septentrio PolaRxS
TRM41249.00
Færøyene
FAR2
1012
Septentrio PolaRxS
LEIAT504GG
Kautokeino
KAU2
1006
Septentrio PolaRxS
TRM559800.00
Honningsvag
HON2
1005
Septentrio PolaRxS
TRM559800.00
Bjørnøya
BJO2
1008
Septentrio PolaRxS
TRM41249.00
Hopen
HOP2
1009
Septentrio PolaRxS
TRM41249.00Table 2
List of NRISN distance differences (km).
Sites
TRO2
VEG2
NYA2
HOF2
FAR2
KAU2
HON2
BJO2
HOP2
TRO2
X
533.5
1053
1578
1433
176
301
540
787
VEG2
—
X
1476
1278
1004
603
822
1018
1288
NYA2
—
—
X
1846
1982
1152
970
526
411
HOF2
—
—
—
X
493.5
1734
1849
1715
1954
FAR2
—
—
—
—
X
1559
1728
1720
1988
KAU2
—
—
—
—
—
X
245
626
837
HON2
—
—
—
—
—
—
X
454
617
BJO2
—
—
—
—
—
—
—
X
280
HOP2
—
—
—
—
—
—
—
—
XTable 3
NRISN ellipsoidal station heights (WGS84).
Site
Height (m)
Tromso
132.4342
Vega
56.4121
Ny-Ålesund
81.2904
Honningsvag
73.8241
Hopen
37.4744
Hofn
82.8036
Færøyene
110.2158
Kautokeino
413.3949
Bjørnøya
53.6580Figure 1
Sites characteristics of Norwegian Ionospheric Scintillation Network.To accurately determine the relative performances of the various algorithms, it is important to test different values for the filter parameters. This also allows us to determine the optimal filter parameters for the computation of the amplitude and phase variances for GNSS signals. In this case, the relevant parameters are the satellite elevation angle, cutoff frequency, and averaging period. The algorithms were tested with a cutoff frequency of0.1
±
0.05 Hz for the signal amplitude and 0.3
±
0.1 Hz for the phase and an averaging interval of 30, 60, or 90 seconds.


## 3. Digital Filters Construction

To study the influence of scintillation on GNSS signals, the standard deviations of the amplitude fluctuationsS
4 and phase fluctuations σ
ϕ must be computed. To achieve this goal, digital filters that satisfy our requirements have to be constructed. For this, we require both a suitable low-pass filter to compute the amplitude scintillation index S
4 and a high-pass filter to compute the phase scintillation index σ
ϕ.Digital Finite Impulse Response (FIR) and Infinite Impulse Response (IIR) filters are chosen and used to construct sixth-order Butterworth, Chebyshev, and Elliptic low-pass filters. Further details on the design and implementation of digital FIR/IIR filters are well described in [38–41]. Once the filter specifications have been determined, the next step is to compute the corresponding filter coefficients in the frequency domain. This is usually done with the computationally efficient Parks-McClellan algorithm. A detailed description of the algorithm is beyond the scope of this paper. The interested reader is referred to [41, Section 7.3]. Since an ideal low-pass filter is impossible to obtain, we often use approximations. In this paper, three different approximations are considered, namely, the Butterworth, Chebyshev, and Elliptic filters.


### 3.1. Cutoff Frequency Determination

A raw GNSS scintillation signal is a nonstationary signal that is affected by the Doppler shift caused by satellite-receiver relative motion (trend), the slowly varying background ionosphere, and scintillation. In order to study the scintillation, trends and slow fluctuations in the raw signal should be removed; that is, any signal components that do not originate from scintillation must be excluded. There exist many methods for carrying out this detrending operation, for example, linear regression by fitting a straight line to the original data or the ratio and difference detrending methods. Detrending is the statistical operation performed to remove the long-term change in the mean function and is equivalent to high-pass filter. That is, the variance at lower frequencies is reduced compared to variance at high frequencies.In general, the data detrending is a preprocessing step to achieve stationarity of the observed signal. The next step is then to analyze the signal statistically, by, for example, computing the first and the second moments of the signal, namely, the phase scintillation indexσ
ϕ and amplitude scintillation index S
4. To perform a proper detrending, the cutoff frequency must be selected properly. The most important part of the scintillation spectrum is around the Fresnel frequency. It also extends to higher frequencies, but with decreasing power per bandwidth. The Fresnel frequency is f
F
=
v
/
2
λ
z [3, 30], where v is the relative velocity between the satellite and ionosphere, λ is the wavelength of the signal (~19 cm), and z denotes the distance between the receiver and the ionosphere (~350 km at 90 degrees elevation). Ideally, the cutoff frequency should be slightly below the Fresnel frequency but still high enough to remove the unwanted signal components. Also, it should be a constant, so that the scintillation values are comparable. For the case of GNSS scintillation data, a Butterworth filter with a constant cutoff frequency of 0.1 Hz is the most common approach, but it has been shown that this value is inappropriate at high latitudes [16, 24].


### 3.2. Butterworth Approximation

The high-order Butterworth approximation [40, p. 264] satisfies our needs, as the filter has a low ripple and minimum transition band. Due to these characteristics, the Butterworth sixth-order filter is the most used filter in computation of scintillation indices S
4 and σ
ϕ. The cutoff frequency f
c is user defined parameter, by default set to 0.1 Hz for low-pass filters and 0.3 Hz for high-pass filters. The derivation of such indices is usually carried out using Algorithm 1, that is, fast convolution using a fast Fourier transform (FFT) and multiplication in the frequency domain. For each Fourier coefficient, the magnitude of the frequency response of nth-order Butterworth filter is given by the following equation:(1)
H
f
=
1
1
+
f
/
f
c
±
2
n
.Note the ± sign in the exponent, which toggles whether the filter exhibits a low-pass or high-pass behavior. Also, the filter gain has been normalized to unity in the region with minimum attenuation, that is, when either f
≪
f
c or f
≫
f
c.


### 3.3. Chebyshev Approximation

The Chebyshev approximation [41, p. 27] has passband ripple that can be dumped and a remarkable transition band. Based on these qualities, its chosen to approximate the ideal low-pass filter with the user specification parameters. However, unlike the Butterworth case, the Chebyshev filter has ripples in the passband and rapid transition and is able to satisfy the user specification with lower order than Butterworth. The phase response is not linear as the Butterworth case.The magnitude of the frequency response of annth-order Chebyshev filter is(2)
H
f
=
1
1
+
ϵ
2
T
n
2
f
/
f
c
,where ϵ determines the ripple magnitude and is less than 1. T
n is a Chebyshev polynomial given by the expression(3)
T
n
x
=
cos
⁡
n
cos
-
1
⁡
x
,
if
x
≤
0
,
cosh
⁡
n
cosh
-
1
⁡
x
,
if
x
>
0
.


### 3.4. Elliptic Approximation

The Elliptic filters [40, p. 275] known asCauer filters yield smaller transition bandwidths than Butterworth or Chebyshev filters for any given order but are equiripple in both the passbands and stopbands. In general, Elliptic filters meet a given set of performance specifications with the lowest order of any filter type.For Elliptic filters, the normalized cutoff frequencyf
c is a number between 0 and 1, where 1 corresponds to half the sampling frequency (Nyquist frequency). The filter has the smallest transition band of the two approaches described so far but is very complicated analytically.The magnitude of the frequency response of annth-order Elliptic filter is(4)
H
f
=
1
1
+
ϵ
2
E
n
2
f
/
f
c
,where ϵ controls the deviation in the passband and E
n is the Jacobian Elliptic function [42, pp. 567–588] of order n.


### 3.5. Implementation

In this study, we have implemented digital filters of types FIR and IIR. The frequency response for all three filter types are given, respectively, by Figures2 and 3. These figures show the difference between the filter types (Butterworth, Elliptic, and Chebyshev).Figure 2
Typical low-pass digital filter characteristics used to compute the amplitude scintillation.Figure 3
Typical high-pass digital filter characteristics used to compute the phase scintillation.


## 3.1. Cutoff Frequency Determination

A raw GNSS scintillation signal is a nonstationary signal that is affected by the Doppler shift caused by satellite-receiver relative motion (trend), the slowly varying background ionosphere, and scintillation. In order to study the scintillation, trends and slow fluctuations in the raw signal should be removed; that is, any signal components that do not originate from scintillation must be excluded. There exist many methods for carrying out this detrending operation, for example, linear regression by fitting a straight line to the original data or the ratio and difference detrending methods. Detrending is the statistical operation performed to remove the long-term change in the mean function and is equivalent to high-pass filter. That is, the variance at lower frequencies is reduced compared to variance at high frequencies.In general, the data detrending is a preprocessing step to achieve stationarity of the observed signal. The next step is then to analyze the signal statistically, by, for example, computing the first and the second moments of the signal, namely, the phase scintillation indexσ
ϕ and amplitude scintillation index S
4. To perform a proper detrending, the cutoff frequency must be selected properly. The most important part of the scintillation spectrum is around the Fresnel frequency. It also extends to higher frequencies, but with decreasing power per bandwidth. The Fresnel frequency is f
F
=
v
/
2
λ
z [3, 30], where v is the relative velocity between the satellite and ionosphere, λ is the wavelength of the signal (~19 cm), and z denotes the distance between the receiver and the ionosphere (~350 km at 90 degrees elevation). Ideally, the cutoff frequency should be slightly below the Fresnel frequency but still high enough to remove the unwanted signal components. Also, it should be a constant, so that the scintillation values are comparable. For the case of GNSS scintillation data, a Butterworth filter with a constant cutoff frequency of 0.1 Hz is the most common approach, but it has been shown that this value is inappropriate at high latitudes [16, 24].


## 3.2. Butterworth Approximation

The high-order Butterworth approximation [40, p. 264] satisfies our needs, as the filter has a low ripple and minimum transition band. Due to these characteristics, the Butterworth sixth-order filter is the most used filter in computation of scintillation indices S
4 and σ
ϕ. The cutoff frequency f
c is user defined parameter, by default set to 0.1 Hz for low-pass filters and 0.3 Hz for high-pass filters. The derivation of such indices is usually carried out using Algorithm 1, that is, fast convolution using a fast Fourier transform (FFT) and multiplication in the frequency domain. For each Fourier coefficient, the magnitude of the frequency response of nth-order Butterworth filter is given by the following equation:(1)
H
f
=
1
1
+
f
/
f
c
±
2
n
.Note the ± sign in the exponent, which toggles whether the filter exhibits a low-pass or high-pass behavior. Also, the filter gain has been normalized to unity in the region with minimum attenuation, that is, when either f
≪
f
c or f
≫
f
c.


## 3.3. Chebyshev Approximation

The Chebyshev approximation [41, p. 27] has passband ripple that can be dumped and a remarkable transition band. Based on these qualities, its chosen to approximate the ideal low-pass filter with the user specification parameters. However, unlike the Butterworth case, the Chebyshev filter has ripples in the passband and rapid transition and is able to satisfy the user specification with lower order than Butterworth. The phase response is not linear as the Butterworth case.The magnitude of the frequency response of annth-order Chebyshev filter is(2)
H
f
=
1
1
+
ϵ
2
T
n
2
f
/
f
c
,where ϵ determines the ripple magnitude and is less than 1. T
n is a Chebyshev polynomial given by the expression(3)
T
n
x
=
cos
⁡
n
cos
-
1
⁡
x
,
if
x
≤
0
,
cosh
⁡
n
cosh
-
1
⁡
x
,
if
x
>
0
.


## 3.4. Elliptic Approximation

The Elliptic filters [40, p. 275] known asCauer filters yield smaller transition bandwidths than Butterworth or Chebyshev filters for any given order but are equiripple in both the passbands and stopbands. In general, Elliptic filters meet a given set of performance specifications with the lowest order of any filter type.For Elliptic filters, the normalized cutoff frequencyf
c is a number between 0 and 1, where 1 corresponds to half the sampling frequency (Nyquist frequency). The filter has the smallest transition band of the two approaches described so far but is very complicated analytically.The magnitude of the frequency response of annth-order Elliptic filter is(4)
H
f
=
1
1
+
ϵ
2
E
n
2
f
/
f
c
,where ϵ controls the deviation in the passband and E
n is the Jacobian Elliptic function [42, pp. 567–588] of order n.


## 3.5. Implementation

In this study, we have implemented digital filters of types FIR and IIR. The frequency response for all three filter types are given, respectively, by Figures2 and 3. These figures show the difference between the filter types (Butterworth, Elliptic, and Chebyshev).Figure 2
Typical low-pass digital filter characteristics used to compute the amplitude scintillation.Figure 3
Typical high-pass digital filter characteristics used to compute the phase scintillation.


## 4. Nonparametric Regression

In contrast to other data modeling methods, for instance, parametric and semiparametric, a nonparametric regression makes no assumptions about the shape of the distribution function. This class of modeling lets the data speak for themselves and they open the way for a new model by their flexibility.Let us now suppose that we have a vector of observationsX
=
(
X
1
,
X
2
,
…
,
X
n
) sampled from an unknown density function f
(
x
) and that our aim is to estimate f
(
x
) and display it graphically. The density estimation function is then defined as(5)
f
^
x
=
1
n
h
∑
i
=
1
n
K
x
i
-
x
h
,where h is thesmoothing parameter and K
(
·
) is the kernel function. This kernel cannot be arbitrarily chosen but needs to satisfy three requirements: it must be nonnegative, L
2 normalized, and symmetric around the origin.There are lots of different kernel functions that satisfy these requirements, for example, the uniform kernel, quartic kernel, Gaussian kernel, and Epanechnikov kernel. The precise choice of kernel function is in practice less important than the choice of smoothing parameterh.


### 4.1. MSE and MISE

We wish to investigate the performance of the kernel density estimation (KDE) at a single point or over the whole real line and find out how close our estimator is to its target. The Mean Square Error (MSE) and Mean Integrated Square Error (MISE) can be used to measure such performance or efficiency of nonparametric methods.The MSE when estimatingf
(
x
) using the estimator f
^
(
x
) at point x is given by(6)
MSE
⁡
f
^
x
=
E
f
^
x
-
f
x
2
=
var
⁡
f
^
x
+
bias
⁡
f
^
x
2
.We see that the MSE can be decomposed into two parts, namely, variance and bias terms. In contrast to parametric models, the bias is ignored. The MISE is defined as(7)
MISE
⁡
f
^
x
=
∫
MSE
⁡
f
^
x
d
x
.


### 4.2. Asymptotic Statistical Properties of KDE

In this section, we will derive asymptotic approximations for the bias and variance of the KDE. These terms are necessary to derive the optimal smoothing parameterh
o
p
t. For these derivations, it will be useful for defining the second moment μ
2
(
K
) and norm K of the KDE: (8)
μ
2
K
=
∫
-
∞
∞
u
2
K
u
d
u
,
K
2
2
=
∫
-
∞
∞
K
u
2
d
u
.


#### 4.2.1. Derivation of Bias Approximation

We will start by deriving an expression for the expectation value of the kernel densityf
^
h
(
x
):(9)
E
f
^
h
x
=
1
n
h
∑
i
=
1
n
E
K
x
-
X
i
h
≈
f
x
+
1
2
h
2
f
′
′
x
μ
2
K
.The approximation ofE
{
f
^
h
(
x
)
} is obtained by using the Taylor expansion up to the second terms in h and letting h
→
0. By the definition of bias, this result yields(10)
bias
⁡
f
^
h
x
=
E
f
^
h
x
-
f
x
≈
1
2
h
2
f
′
′
x
μ
2
K
.Some facts about the bias off
^
h
(
x
) are as follows:(1)
bias
⁡
f
^
h
x
~
o
(
h
2
). Larger values of h will characterize the process as oversmoothing.
(2)
The sign and the direction are decided byf
′
′
(
x
).
(3)
If we include more terms in the Taylor series expansion, we will reduce the error term.


#### 4.2.2. Derivation of Variance Approximation

The variance of stochastic variableX is defined as var
⁡
(
X
)
=
E
(
X
2
)
-
E
(
X
)
2 and applying the Taylor expansion of f
(
x
+
s
h
) and ∫
-
∞
∞
‍
s
K
(
s
)
2
d
s
=
0 (from the symmetry of the kernel):(11)
var
⁡
f
^
h
x
=
1
n
2
var
⁡
1
h
∑
i
=
1
n
K
x
-
X
i
h
=
1
n
h
K
2
2
f
x
+
o
1
n
h
.Some facts about the variance off
^
h
(
x
) are as follows:(1)
var
⁡
f
^
h
x
~
o
(
h
-
1
). Small values of h will characterize the process as undersmoothing.
(2)
n
h gives the number of the observations inside the processing window.
(3)
The variance increases with the magnitude off
(
x
).


#### 4.2.3. Bandwidth Selection of KDE

As mentioned earlier, the bandwidth selection is more important than the choice of the kernel. This section is devoted to computing the optimal smoothing parameterh
o
p
t. We start by evaluating the expressions given by (6) and (7):(12)
MSE
⁡
f
^
h
x
=
E
f
^
h
x
-
f
x
2
≈
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
f
x
,
(13)
MISE
⁡
f
^
h
x
=
E
∫
-
∞
∞
f
^
h
x
-
f
x
2
d
x
≈
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
+
o
n
h
-
1
+
o
h
4for h
→
0 and (
n
h
)
→
∞, and the Asymptotic Mean Integrated Square Error (AMISE) reads(14)
AMISE
⁡
f
^
h
x
=
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
.The exact computation ofM
I
S
E is given in Wand and Jones [43, p. 24].The optimal smoothing parameterh
o
p
t is found by minimizing the expression given by (14):(15)
∂
AMISE
⁡
f
^
h
x
∂
h
=
d
d
h
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
;set (15) to be equal to 0, and solve; we get(16)
h
M
I
S
E
~
K
2
2
f
′
′
x
μ
2
K
2
n
1
/
5
~
o
n
-
1
/
5
.Inserting this expression of h
M
I
S
E into (13), we get(17)
MISE
⁡
f
^
h
x
~
5
4
f
′
′
x
K
4
μ
2
K
2
1
/
5
n
-
4
/
5
~
o
n
-
4
/
5
.Several smoothing parameter selection procedures exist, for instance, plug-in methods (Section4.2.3), focused information criteria (FIC) [44, p. 145], cross-validation, penalizing functions, and Akaike Information Criteria (AIC). The interested reader is referred to [45, Chap. 5].


### 4.3. Kernel Methods for Nonparametric Regression

In this section, we try to put all pieces together to construct a class of kernel-type regression estimators known aslocal polynomial kernel estimators. The main idea is to estimate the regression function at a particular point x by locally fitting a pth degree polynomial to the data by employing the weighted least square techniques. The weights are chosen to the height of the kernel function centered about the point x.The steps needed to derive the expression of the local polynomial kernel estimator are as follows.(i)
Local polynomial definition: Letp be the order of the polynomial we fit at point x to estimate y
^
=
f
^
(
x
):(18)
P
x
=
β
0
+
β
1
X
i
-
x
+
⋯
+
β
p
X
i
-
x
p
.
(ii)
Weights definition: the weights are given by the kernel density function(19)
K
h
X
i
-
x
=
1
h
K
X
i
-
x
h
.
(iii)
Weighted least square: the value of the estimate at a pointx is β
0, where β
^ minimize the expression:(20)
S
β
=
∑
K
h
X
i
-
x
Y
i
-
P
x
2
.
The weighted least square solution in matrix form readβ
^
=
(
X
T
W
X
)
-
1
X
T
W
Y, where Y is (
n
×
1
) vector of responses. The design matrix X read(21)
X
=
1
X
1
-
x
⋯
X
1
-
x
p
⋮
⋮
⋱
⋮
1
X
n
-
x
⋯
X
n
-
x
p
.
The weight matrixW is (
n
×
n
) diagonal matrix of weights w
i
=
K
h
(
X
i
-
x
):(22)
W
=
w
1
0
⋯
0
0
w
2
⋮
0
⋮
⋮
⋱
⋮
0
0
⋯
w
n
.
The value of the estimate at a pointx is the intercept coefficient β
0 of the local fit of P
(
x
):(23)
y
^
=
e
T
X
T
W
X
-
1
X
T
W
Y
,
wheree
T
=
1
0
⋯
0 is a vector of dimensions (
(
p
+
1
)
×
1
).
Some important remarks are as follows.
(1)
Forp
=
0, we fit a constant function locally and this estimator is known as the Nadaraya-Watson kernel estimator:(24)
f
^
N
W
x
=
∑
i
=
1
n
K
h
X
i
-
x
Y
i
∑
i
=
1
n
K
h
X
i
-
x
.
(2)
Forp
=
1, the estimator function corresponds to the Priestley-Chao kernel estimator(25)
f
^
P
C
x
=
1
h
∑
i
=
1
n
x
i
-
x
i
-
1
K
x
-
x
i
h
y
i
.
(3)
Final product: the key is the optimal smoothing parameterh
o
p
t. To compute the phase scintillation index σ
ϕ, the steps are as follows:
(i)
Find the smoothing parameterh
o
p
t. In this study, we have employed the AIC and the bias Corrected AIC (AICC) for smoothing parameter selection [46, Chap. 3]. This operation is repeated for each computation of σ
ϕ.
(ii)
Compute the weights based on the height of the chosen kernel.
(iii)
Locally fit apth degree polynomial to the data by using the weighted least square technique.


### 4.4. Analytical Analysis of KDE

We will now derive an analytical relationship between the conventional filters and the new approach based on kernel smoothing. With the conventional approach, we apply a Fourier transform to the functiong
(
x
), multiply it by the filter transfer function H
(
f
) in the frequency domain, and then apply the inverse Fourier transform to transform the filtered function back to its original domain. According to the convolution theorem of Fourier analysis, this can equivalently be written as a convolution of the function g
(
x
) with the inverse Fourier transform h
(
x
)
=
F
-
1
{
H
(
f
)
} of the filter transfer function:(26)
F
-
1
F
g
x
H
f
=
g
x
∗
h
x
.For simplicity, we will limit the analysis to a first-order low-pass filter of the Butterworth type, as obtained by setting n
=
1 and choosing a positive exponent in (1):(27)
H
f
=
1
1
+
f
/
f
c
2
. Taking the inverse Fourier transform of this, we get(28)
h
x
=
f
c
π
K
0
f
c
x
,where K
0
(
·
) is the modified Bessel function of the second kind. Substituting back into (26), we get(29)
F
-
1
F
g
x
H
f
=
f
c
π
g
x
∗
K
0
f
c
x
. Writing the convolution explicitly out, we find that the first-order Butterworth filter is equivalent to a kernel smoothing using a Bessel function as the kernel:(30)
F
-
1
F
g
x
H
f
=
f
c
π
∫
g
x
K
0
f
c
x
-
x
′
d
x
′
.From this result, we see that the cutoff frequency f
c of the Butterworth filter is directly equivalent to the smoothing parameter h in KDE equation (5). The modified Bessel function of the second kind can be written in integral form as(31)
K
0
z
=
∫
0
∞
cos
⁡
z
t
1
+
t
2
d
t
. It is also worth noting that(32)
1
π
∫
-
∞
∞
K
0
z
d
z
=
1
, which means that the Bessel kernel is a proper density function, but with a pole at the origin. We have established a direct analytical link between the conventional Butterworth filter approach and the new kernel smoothing approach. In this light, the nonparametric filter approach can be considered as a generalization of the conventional filters to other kernels, with the conventional filters as a special case.


#### 4.4.1. Higher Order Approximation

For higher ordern, the computation of the Fourier transforms of (1), (2), and (4) is very complex analytically. For future work, the numerical methods and Monte Carlo simulation may be considered as suitable techniques for analyzing the relationship between the conventional filters and the new approach based on nonparametric local regression.


#### 4.4.2. Scintillation Cutoff Frequency Determination

We have established the link between the optimal smoothing parameter1
/
h
o
p
t
=
B
w and the cutoff frequency f
c of associated filter equations (30) and (5). As a result, we can define the filter link function implicitly by Ψ
=
F
(
f
c
,
B
w
)
=
0. The relation between f
c and B
w is given by the following equation:(33)
f
c
=
κ
B
w
or
B
w
=
κ
f
c
,where κ is a proportionality constant.


## 4.1. MSE and MISE

We wish to investigate the performance of the kernel density estimation (KDE) at a single point or over the whole real line and find out how close our estimator is to its target. The Mean Square Error (MSE) and Mean Integrated Square Error (MISE) can be used to measure such performance or efficiency of nonparametric methods.The MSE when estimatingf
(
x
) using the estimator f
^
(
x
) at point x is given by(6)
MSE
⁡
f
^
x
=
E
f
^
x
-
f
x
2
=
var
⁡
f
^
x
+
bias
⁡
f
^
x
2
.We see that the MSE can be decomposed into two parts, namely, variance and bias terms. In contrast to parametric models, the bias is ignored. The MISE is defined as(7)
MISE
⁡
f
^
x
=
∫
MSE
⁡
f
^
x
d
x
.


## 4.2. Asymptotic Statistical Properties of KDE

In this section, we will derive asymptotic approximations for the bias and variance of the KDE. These terms are necessary to derive the optimal smoothing parameterh
o
p
t. For these derivations, it will be useful for defining the second moment μ
2
(
K
) and norm K of the KDE: (8)
μ
2
K
=
∫
-
∞
∞
u
2
K
u
d
u
,
K
2
2
=
∫
-
∞
∞
K
u
2
d
u
.


### 4.2.1. Derivation of Bias Approximation

We will start by deriving an expression for the expectation value of the kernel densityf
^
h
(
x
):(9)
E
f
^
h
x
=
1
n
h
∑
i
=
1
n
E
K
x
-
X
i
h
≈
f
x
+
1
2
h
2
f
′
′
x
μ
2
K
.The approximation ofE
{
f
^
h
(
x
)
} is obtained by using the Taylor expansion up to the second terms in h and letting h
→
0. By the definition of bias, this result yields(10)
bias
⁡
f
^
h
x
=
E
f
^
h
x
-
f
x
≈
1
2
h
2
f
′
′
x
μ
2
K
.Some facts about the bias off
^
h
(
x
) are as follows:(1)
bias
⁡
f
^
h
x
~
o
(
h
2
). Larger values of h will characterize the process as oversmoothing.
(2)
The sign and the direction are decided byf
′
′
(
x
).
(3)
If we include more terms in the Taylor series expansion, we will reduce the error term.


### 4.2.2. Derivation of Variance Approximation

The variance of stochastic variableX is defined as var
⁡
(
X
)
=
E
(
X
2
)
-
E
(
X
)
2 and applying the Taylor expansion of f
(
x
+
s
h
) and ∫
-
∞
∞
‍
s
K
(
s
)
2
d
s
=
0 (from the symmetry of the kernel):(11)
var
⁡
f
^
h
x
=
1
n
2
var
⁡
1
h
∑
i
=
1
n
K
x
-
X
i
h
=
1
n
h
K
2
2
f
x
+
o
1
n
h
.Some facts about the variance off
^
h
(
x
) are as follows:(1)
var
⁡
f
^
h
x
~
o
(
h
-
1
). Small values of h will characterize the process as undersmoothing.
(2)
n
h gives the number of the observations inside the processing window.
(3)
The variance increases with the magnitude off
(
x
).


### 4.2.3. Bandwidth Selection of KDE

As mentioned earlier, the bandwidth selection is more important than the choice of the kernel. This section is devoted to computing the optimal smoothing parameterh
o
p
t. We start by evaluating the expressions given by (6) and (7):(12)
MSE
⁡
f
^
h
x
=
E
f
^
h
x
-
f
x
2
≈
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
f
x
,
(13)
MISE
⁡
f
^
h
x
=
E
∫
-
∞
∞
f
^
h
x
-
f
x
2
d
x
≈
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
+
o
n
h
-
1
+
o
h
4for h
→
0 and (
n
h
)
→
∞, and the Asymptotic Mean Integrated Square Error (AMISE) reads(14)
AMISE
⁡
f
^
h
x
=
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
.The exact computation ofM
I
S
E is given in Wand and Jones [43, p. 24].The optimal smoothing parameterh
o
p
t is found by minimizing the expression given by (14):(15)
∂
AMISE
⁡
f
^
h
x
∂
h
=
d
d
h
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
;set (15) to be equal to 0, and solve; we get(16)
h
M
I
S
E
~
K
2
2
f
′
′
x
μ
2
K
2
n
1
/
5
~
o
n
-
1
/
5
.Inserting this expression of h
M
I
S
E into (13), we get(17)
MISE
⁡
f
^
h
x
~
5
4
f
′
′
x
K
4
μ
2
K
2
1
/
5
n
-
4
/
5
~
o
n
-
4
/
5
.Several smoothing parameter selection procedures exist, for instance, plug-in methods (Section4.2.3), focused information criteria (FIC) [44, p. 145], cross-validation, penalizing functions, and Akaike Information Criteria (AIC). The interested reader is referred to [45, Chap. 5].


## 4.2.1. Derivation of Bias Approximation

We will start by deriving an expression for the expectation value of the kernel densityf
^
h
(
x
):(9)
E
f
^
h
x
=
1
n
h
∑
i
=
1
n
E
K
x
-
X
i
h
≈
f
x
+
1
2
h
2
f
′
′
x
μ
2
K
.The approximation ofE
{
f
^
h
(
x
)
} is obtained by using the Taylor expansion up to the second terms in h and letting h
→
0. By the definition of bias, this result yields(10)
bias
⁡
f
^
h
x
=
E
f
^
h
x
-
f
x
≈
1
2
h
2
f
′
′
x
μ
2
K
.Some facts about the bias off
^
h
(
x
) are as follows:(1)
bias
⁡
f
^
h
x
~
o
(
h
2
). Larger values of h will characterize the process as oversmoothing.
(2)
The sign and the direction are decided byf
′
′
(
x
).
(3)
If we include more terms in the Taylor series expansion, we will reduce the error term.


## 4.2.2. Derivation of Variance Approximation

The variance of stochastic variableX is defined as var
⁡
(
X
)
=
E
(
X
2
)
-
E
(
X
)
2 and applying the Taylor expansion of f
(
x
+
s
h
) and ∫
-
∞
∞
‍
s
K
(
s
)
2
d
s
=
0 (from the symmetry of the kernel):(11)
var
⁡
f
^
h
x
=
1
n
2
var
⁡
1
h
∑
i
=
1
n
K
x
-
X
i
h
=
1
n
h
K
2
2
f
x
+
o
1
n
h
.Some facts about the variance off
^
h
(
x
) are as follows:(1)
var
⁡
f
^
h
x
~
o
(
h
-
1
). Small values of h will characterize the process as undersmoothing.
(2)
n
h gives the number of the observations inside the processing window.
(3)
The variance increases with the magnitude off
(
x
).


## 4.2.3. Bandwidth Selection of KDE

As mentioned earlier, the bandwidth selection is more important than the choice of the kernel. This section is devoted to computing the optimal smoothing parameterh
o
p
t. We start by evaluating the expressions given by (6) and (7):(12)
MSE
⁡
f
^
h
x
=
E
f
^
h
x
-
f
x
2
≈
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
f
x
,
(13)
MISE
⁡
f
^
h
x
=
E
∫
-
∞
∞
f
^
h
x
-
f
x
2
d
x
≈
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
+
o
n
h
-
1
+
o
h
4for h
→
0 and (
n
h
)
→
∞, and the Asymptotic Mean Integrated Square Error (AMISE) reads(14)
AMISE
⁡
f
^
h
x
=
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
.The exact computation ofM
I
S
E is given in Wand and Jones [43, p. 24].The optimal smoothing parameterh
o
p
t is found by minimizing the expression given by (14):(15)
∂
AMISE
⁡
f
^
h
x
∂
h
=
d
d
h
h
4
4
f
′
′
x
μ
2
K
2
+
1
n
h
K
2
2
;set (15) to be equal to 0, and solve; we get(16)
h
M
I
S
E
~
K
2
2
f
′
′
x
μ
2
K
2
n
1
/
5
~
o
n
-
1
/
5
.Inserting this expression of h
M
I
S
E into (13), we get(17)
MISE
⁡
f
^
h
x
~
5
4
f
′
′
x
K
4
μ
2
K
2
1
/
5
n
-
4
/
5
~
o
n
-
4
/
5
.Several smoothing parameter selection procedures exist, for instance, plug-in methods (Section4.2.3), focused information criteria (FIC) [44, p. 145], cross-validation, penalizing functions, and Akaike Information Criteria (AIC). The interested reader is referred to [45, Chap. 5].


## 4.3. Kernel Methods for Nonparametric Regression

In this section, we try to put all pieces together to construct a class of kernel-type regression estimators known aslocal polynomial kernel estimators. The main idea is to estimate the regression function at a particular point x by locally fitting a pth degree polynomial to the data by employing the weighted least square techniques. The weights are chosen to the height of the kernel function centered about the point x.The steps needed to derive the expression of the local polynomial kernel estimator are as follows.(i)
Local polynomial definition: Letp be the order of the polynomial we fit at point x to estimate y
^
=
f
^
(
x
):(18)
P
x
=
β
0
+
β
1
X
i
-
x
+
⋯
+
β
p
X
i
-
x
p
.
(ii)
Weights definition: the weights are given by the kernel density function(19)
K
h
X
i
-
x
=
1
h
K
X
i
-
x
h
.
(iii)
Weighted least square: the value of the estimate at a pointx is β
0, where β
^ minimize the expression:(20)
S
β
=
∑
K
h
X
i
-
x
Y
i
-
P
x
2
.
The weighted least square solution in matrix form readβ
^
=
(
X
T
W
X
)
-
1
X
T
W
Y, where Y is (
n
×
1
) vector of responses. The design matrix X read(21)
X
=
1
X
1
-
x
⋯
X
1
-
x
p
⋮
⋮
⋱
⋮
1
X
n
-
x
⋯
X
n
-
x
p
.
The weight matrixW is (
n
×
n
) diagonal matrix of weights w
i
=
K
h
(
X
i
-
x
):(22)
W
=
w
1
0
⋯
0
0
w
2
⋮
0
⋮
⋮
⋱
⋮
0
0
⋯
w
n
.
The value of the estimate at a pointx is the intercept coefficient β
0 of the local fit of P
(
x
):(23)
y
^
=
e
T
X
T
W
X
-
1
X
T
W
Y
,
wheree
T
=
1
0
⋯
0 is a vector of dimensions (
(
p
+
1
)
×
1
).
Some important remarks are as follows.
(1)
Forp
=
0, we fit a constant function locally and this estimator is known as the Nadaraya-Watson kernel estimator:(24)
f
^
N
W
x
=
∑
i
=
1
n
K
h
X
i
-
x
Y
i
∑
i
=
1
n
K
h
X
i
-
x
.
(2)
Forp
=
1, the estimator function corresponds to the Priestley-Chao kernel estimator(25)
f
^
P
C
x
=
1
h
∑
i
=
1
n
x
i
-
x
i
-
1
K
x
-
x
i
h
y
i
.
(3)
Final product: the key is the optimal smoothing parameterh
o
p
t. To compute the phase scintillation index σ
ϕ, the steps are as follows:
(i)
Find the smoothing parameterh
o
p
t. In this study, we have employed the AIC and the bias Corrected AIC (AICC) for smoothing parameter selection [46, Chap. 3]. This operation is repeated for each computation of σ
ϕ.
(ii)
Compute the weights based on the height of the chosen kernel.
(iii)
Locally fit apth degree polynomial to the data by using the weighted least square technique.


## 4.4. Analytical Analysis of KDE

We will now derive an analytical relationship between the conventional filters and the new approach based on kernel smoothing. With the conventional approach, we apply a Fourier transform to the functiong
(
x
), multiply it by the filter transfer function H
(
f
) in the frequency domain, and then apply the inverse Fourier transform to transform the filtered function back to its original domain. According to the convolution theorem of Fourier analysis, this can equivalently be written as a convolution of the function g
(
x
) with the inverse Fourier transform h
(
x
)
=
F
-
1
{
H
(
f
)
} of the filter transfer function:(26)
F
-
1
F
g
x
H
f
=
g
x
∗
h
x
.For simplicity, we will limit the analysis to a first-order low-pass filter of the Butterworth type, as obtained by setting n
=
1 and choosing a positive exponent in (1):(27)
H
f
=
1
1
+
f
/
f
c
2
. Taking the inverse Fourier transform of this, we get(28)
h
x
=
f
c
π
K
0
f
c
x
,where K
0
(
·
) is the modified Bessel function of the second kind. Substituting back into (26), we get(29)
F
-
1
F
g
x
H
f
=
f
c
π
g
x
∗
K
0
f
c
x
. Writing the convolution explicitly out, we find that the first-order Butterworth filter is equivalent to a kernel smoothing using a Bessel function as the kernel:(30)
F
-
1
F
g
x
H
f
=
f
c
π
∫
g
x
K
0
f
c
x
-
x
′
d
x
′
.From this result, we see that the cutoff frequency f
c of the Butterworth filter is directly equivalent to the smoothing parameter h in KDE equation (5). The modified Bessel function of the second kind can be written in integral form as(31)
K
0
z
=
∫
0
∞
cos
⁡
z
t
1
+
t
2
d
t
. It is also worth noting that(32)
1
π
∫
-
∞
∞
K
0
z
d
z
=
1
, which means that the Bessel kernel is a proper density function, but with a pole at the origin. We have established a direct analytical link between the conventional Butterworth filter approach and the new kernel smoothing approach. In this light, the nonparametric filter approach can be considered as a generalization of the conventional filters to other kernels, with the conventional filters as a special case.


### 4.4.1. Higher Order Approximation

For higher ordern, the computation of the Fourier transforms of (1), (2), and (4) is very complex analytically. For future work, the numerical methods and Monte Carlo simulation may be considered as suitable techniques for analyzing the relationship between the conventional filters and the new approach based on nonparametric local regression.


### 4.4.2. Scintillation Cutoff Frequency Determination

We have established the link between the optimal smoothing parameter1
/
h
o
p
t
=
B
w and the cutoff frequency f
c of associated filter equations (30) and (5). As a result, we can define the filter link function implicitly by Ψ
=
F
(
f
c
,
B
w
)
=
0. The relation between f
c and B
w is given by the following equation:(33)
f
c
=
κ
B
w
or
B
w
=
κ
f
c
,where κ is a proportionality constant.


## 4.4.1. Higher Order Approximation

For higher ordern, the computation of the Fourier transforms of (1), (2), and (4) is very complex analytically. For future work, the numerical methods and Monte Carlo simulation may be considered as suitable techniques for analyzing the relationship between the conventional filters and the new approach based on nonparametric local regression.


## 4.4.2. Scintillation Cutoff Frequency Determination

We have established the link between the optimal smoothing parameter1
/
h
o
p
t
=
B
w and the cutoff frequency f
c of associated filter equations (30) and (5). As a result, we can define the filter link function implicitly by Ψ
=
F
(
f
c
,
B
w
)
=
0. The relation between f
c and B
w is given by the following equation:(33)
f
c
=
κ
B
w
or
B
w
=
κ
f
c
,where κ is a proportionality constant.


## 5. Scintillation Indices Computation

This section presents the methods used to compute the amplitude scintillation indexS
4 and phase scintillation index σ
ϕ, that is, to suppress all contributions but the scintillation effects in the signal variance. Following the mathematical formulation and notation suggested by Kaplan et al. [47, chap. 6.4], a simplified model for a received signal in the absence of scintillation is(34)
r
t
=
2
P
s
t
cos
⁡
ω
t
+
ϕ
+
n
t
,where P is the received signal power, r
(
t
) is the observed signal, s
(
t
) is the normalized transmitted signal, n
(
t
) is the signal noise, t is the time, ω is the carrier frequency, and ϕ is a phase offset. Scintillation causes a perturbation to both the received signal amplitude and the phase. Thus, in the presence of scintillation, the model may be extended as(35)
r
t
=
2
P
δ
P
s
t
cos
⁡
ω
t
+
ϕ
+
δ
ϕ
+
n
t
,where δ
P is a positive number that parametrizes signal attenuation, while δ
ϕ represents fluctuations in the phase offset.Figure4 shows the steps necessary to compute the indices discussed above.Figure 4
Building blocks of scintillation indices computation.The first step in the processing chain is the data cleaning, which includes detecting and pruning outliers and handling discontinuities caused by cycle-slips and missing observations in the data set. Next, we have to detrend the data set. This means that we construct a time series of the received signal, estimate the trend, and remove this from the signal. This will remove disturbances due to, for example, the Doppler shift from the relative motion of satellite and receiver. Next, we pass the signal through either a high-pass or low-pass filter, in order to isolate the scintillation components in the signal. Finally, the last step is the computation of the scintillation indices, which is the interesting part.


### 5.1. Filtering Techniques

Using digital filters to extract specific signal components is an important technique in many practical applications. For instance, we need the low-frequency components of the signal strength to calculate the amplitude scintillationS
4 and the high-frequency components of the signal phase to calculate σ
ϕ.There are a variety of methods that can be used to design the digital filters. One commonly used method is the bilinear transform. The method uses the analog filter approximation functions that have been developed and translate the filter coefficients in a way that will make them usable for discrete systems, the IIR filter. That is, the output of the filter will depend on the previous values of the output as well as on the past and the current values of the input. Such a process is known as the autoregressive moving average (ARMA) process.Other methods use frequency response of the desired filter to directly determine the digital filter coefficients. These are known as Finite Impulse Response (FIR) filters and are implemented as moving average filters. Once the FIR filter coefficients have been determined, the output signal is simply the convolution of the input signal with the filter coefficients.One efficient way to implement this approach computationally is that the convolution is performed through multiplication in the frequency domain. This is done by first using a fast Fourier transform (FFT) and then multiplying the filter transfer function with the input signal and performing the inverse transform. This is described by Algorithm1, which assumes that the filter coefficients have already been determined.Algorithm 1 (filtering with FFT convolution).
(
1) Apply the Fourier transform to the filter coefficients h
(
n
):(36)
H
z
=
F
F
T
h
n
.
(2) Apply the Fourier transform to the input signal x
(
n
):(37)
X
z
=
F
F
T
x
n
.
(3) Multiply the complex sequences elementwise:(38)
Y
z
=
H
z
X
z
.
(4) Apply the inverse Fourier transform to transform the filtered function back to its original domain:(39)
y
n
=
F
F
T
-
1
Y
z
.
(5) Normalization: The IFFT coefficients are the complex conjugates of the Discrete Fast Transform (DFT) coefficients. The output signal has to be scaled down by the length of the signal.


### 5.2. Computation of the Amplitude Scintillation IndexS
4

The amplitude scintillation indexS
4 is defined as the standard deviation of the signal power normalized to the mean value over the interval of interest, usually a 1-minute interval.In order to compute this index, the necessary steps are as follows:(I)
Computation of the total amplitude indexS
4
T. The index S
4
T is defined by the expression(40)
S
4
T
=
S
I
2
-
S
I
2
S
I
2
,
whereS
I is the detrended power or satellites signal intensity.
(II)
Correction ofS
4 by removing the ambient noise: S
4
T can have large values due to the ambient noise, and this must be removed. The removal process is done by estimating the average carrier-to-noise ratio C
/
N
0 density over the entire evaluation interval and using the estimate to compute S
4 due the ambient noise. If the signal-to-noise (S
/
N) density is known, the predicted S
4
N due to the ambient noise [48] is given by the expression(41)
S
4
N
=
100
S
/
N
1
+
500
19
S
/
N
.
ReplacingS
/
N with the 60-second estimate S
^
/
N
^, we obtain an estimate of S
4 due to the noise S
4
N. Subtracting the square of this value from the square of (40) yields the revised estimate of S
4:(42)
S
4
=
S
I
2
-
S
I
2
S
I
2
-
100
S
^
/
N
^
1
+
500
19
S
^
/
N
^
.
(III)
Detrending the signal to achieve stationarity: Detrending the raw signal intensity is accomplished by dividing the original time series by the filtered one. The effect of this is the removal of the low-frequency variance. The ratio method is attractive because it is dimensionless; however, it cannot be used if the data contains negative values, and it can also become problematic if the fitted line crosses zero.
Detrending by the difference method produces stationarity of the signal under study. The first difference∇
x
t
=
x
t
-
x
t
-
1 eliminates a linear trend, while the second difference ∇
2
x
t
=
x
t
-
2
x
t
-
1
+
x
t
-
2 can eliminate a quadratic trend and so on. In this study, both methods are used and they produce almost the same result; the first difference method has a little higher range than the ratio method.
In order to classify the scintillation activity indicated byS
4, four categories [49] are defined: S
4
≤
0.25 is quiet, S
4
∈
(
0.25,0.50
] is moderate, S
4
∈
(
0.50,1.00
] is disturbed, and S
4
>
1.00 is severe.
(IV)
Practical considerations: A full understanding of the difference betweenS
/
N
0 and C
/
N
0 is necessary and useful for the users of the GNSS receivers. S
/
N
0is usually expressed in dB; it refers to the ratio of the signal power and noise power in the given bandwidth while C
/
N
0 is expressed in dB-Hz; it refers to carrier power and the noise power per unit bandwidth.
The signal strength indicator measured by the receiver isC
/
N
0, the carrier-to-noise density measured in dB-Hz. In the computation of S
4
T, the signal-to-noise density S
/
N
0 is used. The relation between these two quantities [50] is given by the expression(43)
S
/
N
0
=
10
C
/
10
N
0
.


### 5.3. Computation of Scintillation Indexσ
ϕ

Computation of the phase scintillationσ
ϕ is straightforward, as the procedure is very similar to the computation of S
4. The main difference exists in data filtering; for σ
ϕ high-pass filter, usually a Butterworth filter with cutoff frequency of 0.3 Hz is used while S
4 is as mentioned in Section 3.2.For typical FIR, IIR, and FFT convolution, the filter transfer function given by (1) is used. The sign in the exponent is negative, which results in the high-pass behavior. Detrending is a necessary step before carrying out the filtering. This is accomplished by any method that is able to produce a stationary time series. In this study, the ratio method is used for data detrending. The next step is filtering. Algorithm 1 described previously is used.For the case of the new approach, the details are given in Section4.3 and are summarized as follows:(I)
Preprocessing step requires handling outliers, missing observations, and cycle-slip detection and repair.
(II)
Choose a kernel density function to compute the weights.
(III)
Model selection with AIC/AICC is to determine the optimal smoothing parameterh
o
p
t. The key is h
o
p
t.
(IV)
Compute the weights based on the height of the chosen kernel centered about the point of interest.
(V)
locally fit apth degree (p
=
0,1,or
2) polynomial to the data by using the weighted least square technique.
(VI)
Computeσ
ϕ by taking the standard deviation of the residuals.


## 5.1. Filtering Techniques

Using digital filters to extract specific signal components is an important technique in many practical applications. For instance, we need the low-frequency components of the signal strength to calculate the amplitude scintillationS
4 and the high-frequency components of the signal phase to calculate σ
ϕ.There are a variety of methods that can be used to design the digital filters. One commonly used method is the bilinear transform. The method uses the analog filter approximation functions that have been developed and translate the filter coefficients in a way that will make them usable for discrete systems, the IIR filter. That is, the output of the filter will depend on the previous values of the output as well as on the past and the current values of the input. Such a process is known as the autoregressive moving average (ARMA) process.Other methods use frequency response of the desired filter to directly determine the digital filter coefficients. These are known as Finite Impulse Response (FIR) filters and are implemented as moving average filters. Once the FIR filter coefficients have been determined, the output signal is simply the convolution of the input signal with the filter coefficients.One efficient way to implement this approach computationally is that the convolution is performed through multiplication in the frequency domain. This is done by first using a fast Fourier transform (FFT) and then multiplying the filter transfer function with the input signal and performing the inverse transform. This is described by Algorithm1, which assumes that the filter coefficients have already been determined.Algorithm 1 (filtering with FFT convolution).
(
1) Apply the Fourier transform to the filter coefficients h
(
n
):(36)
H
z
=
F
F
T
h
n
.
(2) Apply the Fourier transform to the input signal x
(
n
):(37)
X
z
=
F
F
T
x
n
.
(3) Multiply the complex sequences elementwise:(38)
Y
z
=
H
z
X
z
.
(4) Apply the inverse Fourier transform to transform the filtered function back to its original domain:(39)
y
n
=
F
F
T
-
1
Y
z
.
(5) Normalization: The IFFT coefficients are the complex conjugates of the Discrete Fast Transform (DFT) coefficients. The output signal has to be scaled down by the length of the signal.


## 5.2. Computation of the Amplitude Scintillation IndexS
4

The amplitude scintillation indexS
4 is defined as the standard deviation of the signal power normalized to the mean value over the interval of interest, usually a 1-minute interval.In order to compute this index, the necessary steps are as follows:(I)
Computation of the total amplitude indexS
4
T. The index S
4
T is defined by the expression(40)
S
4
T
=
S
I
2
-
S
I
2
S
I
2
,
whereS
I is the detrended power or satellites signal intensity.
(II)
Correction ofS
4 by removing the ambient noise: S
4
T can have large values due to the ambient noise, and this must be removed. The removal process is done by estimating the average carrier-to-noise ratio C
/
N
0 density over the entire evaluation interval and using the estimate to compute S
4 due the ambient noise. If the signal-to-noise (S
/
N) density is known, the predicted S
4
N due to the ambient noise [48] is given by the expression(41)
S
4
N
=
100
S
/
N
1
+
500
19
S
/
N
.
ReplacingS
/
N with the 60-second estimate S
^
/
N
^, we obtain an estimate of S
4 due to the noise S
4
N. Subtracting the square of this value from the square of (40) yields the revised estimate of S
4:(42)
S
4
=
S
I
2
-
S
I
2
S
I
2
-
100
S
^
/
N
^
1
+
500
19
S
^
/
N
^
.
(III)
Detrending the signal to achieve stationarity: Detrending the raw signal intensity is accomplished by dividing the original time series by the filtered one. The effect of this is the removal of the low-frequency variance. The ratio method is attractive because it is dimensionless; however, it cannot be used if the data contains negative values, and it can also become problematic if the fitted line crosses zero.
Detrending by the difference method produces stationarity of the signal under study. The first difference∇
x
t
=
x
t
-
x
t
-
1 eliminates a linear trend, while the second difference ∇
2
x
t
=
x
t
-
2
x
t
-
1
+
x
t
-
2 can eliminate a quadratic trend and so on. In this study, both methods are used and they produce almost the same result; the first difference method has a little higher range than the ratio method.
In order to classify the scintillation activity indicated byS
4, four categories [49] are defined: S
4
≤
0.25 is quiet, S
4
∈
(
0.25,0.50
] is moderate, S
4
∈
(
0.50,1.00
] is disturbed, and S
4
>
1.00 is severe.
(IV)
Practical considerations: A full understanding of the difference betweenS
/
N
0 and C
/
N
0 is necessary and useful for the users of the GNSS receivers. S
/
N
0is usually expressed in dB; it refers to the ratio of the signal power and noise power in the given bandwidth while C
/
N
0 is expressed in dB-Hz; it refers to carrier power and the noise power per unit bandwidth.
The signal strength indicator measured by the receiver isC
/
N
0, the carrier-to-noise density measured in dB-Hz. In the computation of S
4
T, the signal-to-noise density S
/
N
0 is used. The relation between these two quantities [50] is given by the expression(43)
S
/
N
0
=
10
C
/
10
N
0
.


## 5.3. Computation of Scintillation Indexσ
ϕ

Computation of the phase scintillationσ
ϕ is straightforward, as the procedure is very similar to the computation of S
4. The main difference exists in data filtering; for σ
ϕ high-pass filter, usually a Butterworth filter with cutoff frequency of 0.3 Hz is used while S
4 is as mentioned in Section 3.2.For typical FIR, IIR, and FFT convolution, the filter transfer function given by (1) is used. The sign in the exponent is negative, which results in the high-pass behavior. Detrending is a necessary step before carrying out the filtering. This is accomplished by any method that is able to produce a stationary time series. In this study, the ratio method is used for data detrending. The next step is filtering. Algorithm 1 described previously is used.For the case of the new approach, the details are given in Section4.3 and are summarized as follows:(I)
Preprocessing step requires handling outliers, missing observations, and cycle-slip detection and repair.
(II)
Choose a kernel density function to compute the weights.
(III)
Model selection with AIC/AICC is to determine the optimal smoothing parameterh
o
p
t. The key is h
o
p
t.
(IV)
Compute the weights based on the height of the chosen kernel centered about the point of interest.
(V)
locally fit apth degree (p
=
0,1,or
2) polynomial to the data by using the weighted least square technique.
(VI)
Computeσ
ϕ by taking the standard deviation of the residuals.


## 6. Statistical Analysis

The main objective of this section is to define the distributions of phase and amplitude scintillation indices (S
4 and σ
ϕ), correlation, and sensitivity analysis of the implemented filtering algorithms. Focus will be directed to the new approach, the nonparametric local regression with smoothing parameter selection. Outliers detection and repair, discontinuities, and missing observations are handled in data preparation.


### 6.1. Statistical Distribution of the Scintillation IndexS
4

The power fluctuations ofδ
P given in (35) are generally modeled as a Nakagami m-distribution. The probability density function (PDF) of this distribution is given by (B.1) with mean value of one and variance 1
/
m.Due to the characteristics of the Nakagami distributions,S
4
=
1
/
m will not exceed 2 [47, p. 296]. Figure 5 shows the computed S
4 index for an arbitrary GPS satellite and Figure 6 shows the distribution of S
4.Figure 5
Amplitude scintillation index,S
4 computed for GPS satellite PRN 23, year 2014 and DOY 50, site Tromso. False scintillation(s) due to low elevation angle(s) is (are) given by shadowed rectangle in the figure.Figure 6
Distribution of the amplitude scintillation indexS
4 computed for GLONASS satellite PRN 13, year 2014 and DOY 50, observed at site Tromso.


### 6.2. Distribution ofσ
ϕ

This subsection is devoted to the computation of phase scintillationσ
ϕ. Figure 7 shows the computed phase scintillation index for GPS satellite. Clearly, the heavy tailed distribution family is the appropriate choice. Here, we found that Frechet distribution equation (B.5) models the empirical data very well.Figure 7
Phase scintillation index,σ
ϕ for GPS satellite PRN 23, year 2014 and DOY 50, site Tromso.The motivations for considering the log-normal as an approximation model of the distribution ofσ
ϕ are as follows: First, the log-normal distribution is a positive real-valued function with a heavy tail that can describe the presence of extreme variability in the data. Secondly, the distribution is simple and practical and provides a very good fit to the data. Finally its variance is scaling.Figure8 shows the distribution of phase scintillation indices for an arbitrary GLONASS satellite, PRN 08. Clearly, the distribution follows the heavy tailed distribution (the Frechet and log-normal ones are appropriate fit).Figure 8
Empirical distribution of phase scintillation for an arbitrary GLONASS satellite, PRN 08. Heavy tailed Frechet distribution is an appropriate fit to phase scintillation.


### 6.3. Correlation Analysis

Often, we like to measure the linear predictability of one signalx
t from another one y
t. Assuming that the variance of both time series is finite, the cross-correlation function [51, p.23] can be used for this purpose. This function is defined as(44)
ρ
x
y
s
,
t
=
γ
x
,
y
s
,
t
γ
x
t
,
t
γ
y
s
,
s
,where γ
x
,
y
(
s
,
t
)
=
E
[
(
x
t
-
μ
x
t
)
(
y
s
-
μ
y
s
)
] is the cross-covariance function between two time series, x
t and y
s, respectively. E
(
·
) is the expectation operator. The Cauchy-Schwartz inequality implies γ
x
,
y
s
,
t
2
≤
γ
x
(
t
,
t
)
γ
y
(
s
,
s
).The phase scintillation indexσ
ϕ is the dominant disturbance at high latitudes and is obtained by high-pass filtering. For this reason, we have chosen two implementations of high-pass filtering algorithms (FFT convolution and IIR) for the correlation analysis. All plots have in common the cutoff frequency f
c
=
0.3Hz and the averaging interval 〈
T
¯
〉
=
1 minute. Three scenarios are presented below, using different threshold values T
h
∈
{
0.1,0.15,0.25
} rad.


#### 6.3.1. Scenario 1:σ
ϕ
≤
0.1 rad

For this scenario, we drop all scintillations indicesσ
ϕ above a threshold of 0.1 rad. Figure 9 shows that the algorithms are highly correlated, and the correlation coefficient ρ varies between 78
% and 83
%. The new approach, that is, the nonparametric regression with Akaike information criterion (AIC) [51, pp. 53-54] model selection, performs very well. The algorithm does not take into consideration any knowledge of physical problems, for instance, multipath, Doppler effects due to satellite motion, and cutoff frequency selection. The key is the smoothing parameter h
o
p
t that determines the trade-off between the variance and the bias terms.Figure 9
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 83.4% forσ
ϕ
≤
0.1 rad.
(a)
(b)
(c)In this scenario, we have dropped all observations above the threshold valueT
h
=
0.1 rad. These observations are classified as extreme and correspond to severe ionospheric activities.


#### 6.3.2. Scenario 2:σ
ϕ
≤
0.15 rad

Figure10 shows that the algorithms are still highly correlated and the correlation coefficient ρ has decreased to values between 77.2
% and 79.3
%. The main reason for this reduction is the inclusion of more extremal events that are classified as severe ionospheric activities. Before the vertical lines indicated in Figure 10, the algorithms are highly correlated. After the lines, decorrelation appears between the two algorithms.Figure 10
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 78% forσ
ϕ
≤
0.15 rad.
(a)
(b)
(c)


#### 6.3.3. Scenario 3:σ
ϕ
≤
0.25 rad

In this scenario, we have included allσ
ϕ values (σ
ϕ
>
0.25 rad are treated as outliers and removed). Figure 11 shows that the correlation coefficient has dropped to values between 69 and 74%. The reason is the inclusion of the extremal events classified as severe.Figure 11
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 74% forσ
ϕ
≤
0.25 rad.
(a)
(b)
(c)The vertical line is used to point out/distinguish two classes; the one to the left is highly correlated while the other is uncorrelated and corresponds to the extreme ionospheric activities and possibly noise in the signal.


### 6.4. Sensitivity Analysis

In order to carry out the sensitivity analysis of the implemented algorithms, a reference algorithm has to be chosen. The most common one used is high-pass filter Butterworth with FFT convolution, order 6, and cutoff frequency of0.3 Hz.Let the functionF
=
f
(
x
,
y
,
z
) represents the phase scintillation index σ
ϕ, depending on parameters (
x
,
y
,
z
), where x and y are the cutoff frequency and the averaging interval parameters, respectively. z
=
g
(
α
1
,
α
2
,
…
,
α
k
) represents the remaining parameters, for instance, cycle-slip detection and repair and handling outliers/glitch. We will ignore the parameter z to facilitate analysis.In order to carry out the sensitivity analysis, keeping all parameters constant and varying one parameter at a time are common. This is done by the following equations:(45)
∂
F
∂
x
y
=
c
o
n
s
t
,
∂
F
∂
y
x
=
c
o
n
s
t
.


#### 6.4.1. Scenario 1: Cutoff Frequency Constant

In this scenario, we kept the cutoff frequencyf
c
=
0.3Hz constant and varied the average interval 〈
T
¯
〉
∈
{
30,60,90
} seconds.Figure12 shows the results. The general shape of the time series is very similar, but some of the modes/spikes are not present in all of them. In general using a longer averaging interval smooths the time series, broadening the spikes, reducing their magnitude, and causing some of the spikes to vanish. The longest interval time series does not appear to offer any significant advantage over the middle interval time series. The 30 sec interval time series contains more fine structure than the 60 sec interval time series. In conclusion, the averaging interval should not exceed 60 seconds.Figure 12
The panels show time series ofσ
ϕ calculated from filtered signal phase, using different settings for the σ
ϕ calculation. The data is from GLONASS satellite 25, year 2014, DOY 50, site Honningsvag. (a), (b), and (c) have been calculated for each 30-second, 60-second, and 90-second interval, respectively.
(a)
(b)
(c)


#### 6.4.2. Scenario 2: Averaging Interval Constant

In this scenario, the averaging interval is constant and equal to60 seconds while varying the cutoff frequency f
c
∈
{
0.2,0.3,0.4
} Hz.Figure13 shows the results. The general shape of the time series is very similar, but some of the modes/spikes are not present in all of them. In general using a higher cutoff frequency results in a reduction of the scintillation index value throughout the time series. At time ≈ 340, there is a large spike which is significantly reduced when changing the cutoff frequency from 0.2 to 0.3 Hz and vanishes for a cutoff frequency of 0.4 Hz. The phase scintillation values are clearly strongly dependent on the cutoff frequency, even for a modest change such as the one seen here. This is a known weakness in the standard techniques of σ
ϕ calculation (e.g., [52]).Figure 13
The panels show time series ofσ
ϕ calculated from filtered signal phase, using different settings for the high-pass filter. The data is from GPS satellite 17, year 2014, DOY 50, site Honningsvag. (a), (b), and (c) have been calculated using a cutoff frequency of 0.2, 0.3, and 0.4 Hz, respectively.
(a)
(b)
(c)Taking into account the first and the second sensitivity analysis, the acceptable range for the computation ofσ
ϕ is f
c
∈
[
0.1,0.3
] and 〈
T
¯
〉
∈
[
30,60
]. These ranges give the maximum variations in the data.


## 6.1. Statistical Distribution of the Scintillation IndexS
4

The power fluctuations ofδ
P given in (35) are generally modeled as a Nakagami m-distribution. The probability density function (PDF) of this distribution is given by (B.1) with mean value of one and variance 1
/
m.Due to the characteristics of the Nakagami distributions,S
4
=
1
/
m will not exceed 2 [47, p. 296]. Figure 5 shows the computed S
4 index for an arbitrary GPS satellite and Figure 6 shows the distribution of S
4.Figure 5
Amplitude scintillation index,S
4 computed for GPS satellite PRN 23, year 2014 and DOY 50, site Tromso. False scintillation(s) due to low elevation angle(s) is (are) given by shadowed rectangle in the figure.Figure 6
Distribution of the amplitude scintillation indexS
4 computed for GLONASS satellite PRN 13, year 2014 and DOY 50, observed at site Tromso.


## 6.2. Distribution ofσ
ϕ

This subsection is devoted to the computation of phase scintillationσ
ϕ. Figure 7 shows the computed phase scintillation index for GPS satellite. Clearly, the heavy tailed distribution family is the appropriate choice. Here, we found that Frechet distribution equation (B.5) models the empirical data very well.Figure 7
Phase scintillation index,σ
ϕ for GPS satellite PRN 23, year 2014 and DOY 50, site Tromso.The motivations for considering the log-normal as an approximation model of the distribution ofσ
ϕ are as follows: First, the log-normal distribution is a positive real-valued function with a heavy tail that can describe the presence of extreme variability in the data. Secondly, the distribution is simple and practical and provides a very good fit to the data. Finally its variance is scaling.Figure8 shows the distribution of phase scintillation indices for an arbitrary GLONASS satellite, PRN 08. Clearly, the distribution follows the heavy tailed distribution (the Frechet and log-normal ones are appropriate fit).Figure 8
Empirical distribution of phase scintillation for an arbitrary GLONASS satellite, PRN 08. Heavy tailed Frechet distribution is an appropriate fit to phase scintillation.


## 6.3. Correlation Analysis

Often, we like to measure the linear predictability of one signalx
t from another one y
t. Assuming that the variance of both time series is finite, the cross-correlation function [51, p.23] can be used for this purpose. This function is defined as(44)
ρ
x
y
s
,
t
=
γ
x
,
y
s
,
t
γ
x
t
,
t
γ
y
s
,
s
,where γ
x
,
y
(
s
,
t
)
=
E
[
(
x
t
-
μ
x
t
)
(
y
s
-
μ
y
s
)
] is the cross-covariance function between two time series, x
t and y
s, respectively. E
(
·
) is the expectation operator. The Cauchy-Schwartz inequality implies γ
x
,
y
s
,
t
2
≤
γ
x
(
t
,
t
)
γ
y
(
s
,
s
).The phase scintillation indexσ
ϕ is the dominant disturbance at high latitudes and is obtained by high-pass filtering. For this reason, we have chosen two implementations of high-pass filtering algorithms (FFT convolution and IIR) for the correlation analysis. All plots have in common the cutoff frequency f
c
=
0.3Hz and the averaging interval 〈
T
¯
〉
=
1 minute. Three scenarios are presented below, using different threshold values T
h
∈
{
0.1,0.15,0.25
} rad.


### 6.3.1. Scenario 1:σ
ϕ
≤
0.1 rad

For this scenario, we drop all scintillations indicesσ
ϕ above a threshold of 0.1 rad. Figure 9 shows that the algorithms are highly correlated, and the correlation coefficient ρ varies between 78
% and 83
%. The new approach, that is, the nonparametric regression with Akaike information criterion (AIC) [51, pp. 53-54] model selection, performs very well. The algorithm does not take into consideration any knowledge of physical problems, for instance, multipath, Doppler effects due to satellite motion, and cutoff frequency selection. The key is the smoothing parameter h
o
p
t that determines the trade-off between the variance and the bias terms.Figure 9
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 83.4% forσ
ϕ
≤
0.1 rad.
(a)
(b)
(c)In this scenario, we have dropped all observations above the threshold valueT
h
=
0.1 rad. These observations are classified as extreme and correspond to severe ionospheric activities.


### 6.3.2. Scenario 2:σ
ϕ
≤
0.15 rad

Figure10 shows that the algorithms are still highly correlated and the correlation coefficient ρ has decreased to values between 77.2
% and 79.3
%. The main reason for this reduction is the inclusion of more extremal events that are classified as severe ionospheric activities. Before the vertical lines indicated in Figure 10, the algorithms are highly correlated. After the lines, decorrelation appears between the two algorithms.Figure 10
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 78% forσ
ϕ
≤
0.15 rad.
(a)
(b)
(c)


### 6.3.3. Scenario 3:σ
ϕ
≤
0.25 rad

In this scenario, we have included allσ
ϕ values (σ
ϕ
>
0.25 rad are treated as outliers and removed). Figure 11 shows that the correlation coefficient has dropped to values between 69 and 74%. The reason is the inclusion of the extremal events classified as severe.Figure 11
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 74% forσ
ϕ
≤
0.25 rad.
(a)
(b)
(c)The vertical line is used to point out/distinguish two classes; the one to the left is highly correlated while the other is uncorrelated and corresponds to the extreme ionospheric activities and possibly noise in the signal.


## 6.3.1. Scenario 1:σ
ϕ
≤
0.1 rad

For this scenario, we drop all scintillations indicesσ
ϕ above a threshold of 0.1 rad. Figure 9 shows that the algorithms are highly correlated, and the correlation coefficient ρ varies between 78
% and 83
%. The new approach, that is, the nonparametric regression with Akaike information criterion (AIC) [51, pp. 53-54] model selection, performs very well. The algorithm does not take into consideration any knowledge of physical problems, for instance, multipath, Doppler effects due to satellite motion, and cutoff frequency selection. The key is the smoothing parameter h
o
p
t that determines the trade-off between the variance and the bias terms.Figure 9
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 83.4% forσ
ϕ
≤
0.1 rad.
(a)
(b)
(c)In this scenario, we have dropped all observations above the threshold valueT
h
=
0.1 rad. These observations are classified as extreme and correspond to severe ionospheric activities.


## 6.3.2. Scenario 2:σ
ϕ
≤
0.15 rad

Figure10 shows that the algorithms are still highly correlated and the correlation coefficient ρ has decreased to values between 77.2
% and 79.3
%. The main reason for this reduction is the inclusion of more extremal events that are classified as severe ionospheric activities. Before the vertical lines indicated in Figure 10, the algorithms are highly correlated. After the lines, decorrelation appears between the two algorithms.Figure 10
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 78% forσ
ϕ
≤
0.15 rad.
(a)
(b)
(c)


## 6.3.3. Scenario 3:σ
ϕ
≤
0.25 rad

In this scenario, we have included allσ
ϕ values (σ
ϕ
>
0.25 rad are treated as outliers and removed). Figure 11 shows that the correlation coefficient has dropped to values between 69 and 74%. The reason is the inclusion of the extremal events classified as severe.Figure 11
(a) shows the correlation between the Butterworth sixth-order filter and the IIR filter. (b) shows the correlation between the Butterworth sixth-order filter and an advanced nonparametric regression with kernel smoothing. Finally, (c) shows the correlation between the Butterworth sixth-order filter and nonparametric regression. The plots show high correlation between filtering algorithms, up to 74% forσ
ϕ
≤
0.25 rad.
(a)
(b)
(c)The vertical line is used to point out/distinguish two classes; the one to the left is highly correlated while the other is uncorrelated and corresponds to the extreme ionospheric activities and possibly noise in the signal.


## 6.4. Sensitivity Analysis

In order to carry out the sensitivity analysis of the implemented algorithms, a reference algorithm has to be chosen. The most common one used is high-pass filter Butterworth with FFT convolution, order 6, and cutoff frequency of0.3 Hz.Let the functionF
=
f
(
x
,
y
,
z
) represents the phase scintillation index σ
ϕ, depending on parameters (
x
,
y
,
z
), where x and y are the cutoff frequency and the averaging interval parameters, respectively. z
=
g
(
α
1
,
α
2
,
…
,
α
k
) represents the remaining parameters, for instance, cycle-slip detection and repair and handling outliers/glitch. We will ignore the parameter z to facilitate analysis.In order to carry out the sensitivity analysis, keeping all parameters constant and varying one parameter at a time are common. This is done by the following equations:(45)
∂
F
∂
x
y
=
c
o
n
s
t
,
∂
F
∂
y
x
=
c
o
n
s
t
.


### 6.4.1. Scenario 1: Cutoff Frequency Constant

In this scenario, we kept the cutoff frequencyf
c
=
0.3Hz constant and varied the average interval 〈
T
¯
〉
∈
{
30,60,90
} seconds.Figure12 shows the results. The general shape of the time series is very similar, but some of the modes/spikes are not present in all of them. In general using a longer averaging interval smooths the time series, broadening the spikes, reducing their magnitude, and causing some of the spikes to vanish. The longest interval time series does not appear to offer any significant advantage over the middle interval time series. The 30 sec interval time series contains more fine structure than the 60 sec interval time series. In conclusion, the averaging interval should not exceed 60 seconds.Figure 12
The panels show time series ofσ
ϕ calculated from filtered signal phase, using different settings for the σ
ϕ calculation. The data is from GLONASS satellite 25, year 2014, DOY 50, site Honningsvag. (a), (b), and (c) have been calculated for each 30-second, 60-second, and 90-second interval, respectively.
(a)
(b)
(c)


### 6.4.2. Scenario 2: Averaging Interval Constant

In this scenario, the averaging interval is constant and equal to60 seconds while varying the cutoff frequency f
c
∈
{
0.2,0.3,0.4
} Hz.Figure13 shows the results. The general shape of the time series is very similar, but some of the modes/spikes are not present in all of them. In general using a higher cutoff frequency results in a reduction of the scintillation index value throughout the time series. At time ≈ 340, there is a large spike which is significantly reduced when changing the cutoff frequency from 0.2 to 0.3 Hz and vanishes for a cutoff frequency of 0.4 Hz. The phase scintillation values are clearly strongly dependent on the cutoff frequency, even for a modest change such as the one seen here. This is a known weakness in the standard techniques of σ
ϕ calculation (e.g., [52]).Figure 13
The panels show time series ofσ
ϕ calculated from filtered signal phase, using different settings for the high-pass filter. The data is from GPS satellite 17, year 2014, DOY 50, site Honningsvag. (a), (b), and (c) have been calculated using a cutoff frequency of 0.2, 0.3, and 0.4 Hz, respectively.
(a)
(b)
(c)Taking into account the first and the second sensitivity analysis, the acceptable range for the computation ofσ
ϕ is f
c
∈
[
0.1,0.3
] and 〈
T
¯
〉
∈
[
30,60
]. These ranges give the maximum variations in the data.


## 6.4.1. Scenario 1: Cutoff Frequency Constant

In this scenario, we kept the cutoff frequencyf
c
=
0.3Hz constant and varied the average interval 〈
T
¯
〉
∈
{
30,60,90
} seconds.Figure12 shows the results. The general shape of the time series is very similar, but some of the modes/spikes are not present in all of them. In general using a longer averaging interval smooths the time series, broadening the spikes, reducing their magnitude, and causing some of the spikes to vanish. The longest interval time series does not appear to offer any significant advantage over the middle interval time series. The 30 sec interval time series contains more fine structure than the 60 sec interval time series. In conclusion, the averaging interval should not exceed 60 seconds.Figure 12
The panels show time series ofσ
ϕ calculated from filtered signal phase, using different settings for the σ
ϕ calculation. The data is from GLONASS satellite 25, year 2014, DOY 50, site Honningsvag. (a), (b), and (c) have been calculated for each 30-second, 60-second, and 90-second interval, respectively.
(a)
(b)
(c)


## 6.4.2. Scenario 2: Averaging Interval Constant

In this scenario, the averaging interval is constant and equal to60 seconds while varying the cutoff frequency f
c
∈
{
0.2,0.3,0.4
} Hz.Figure13 shows the results. The general shape of the time series is very similar, but some of the modes/spikes are not present in all of them. In general using a higher cutoff frequency results in a reduction of the scintillation index value throughout the time series. At time ≈ 340, there is a large spike which is significantly reduced when changing the cutoff frequency from 0.2 to 0.3 Hz and vanishes for a cutoff frequency of 0.4 Hz. The phase scintillation values are clearly strongly dependent on the cutoff frequency, even for a modest change such as the one seen here. This is a known weakness in the standard techniques of σ
ϕ calculation (e.g., [52]).Figure 13
The panels show time series ofσ
ϕ calculated from filtered signal phase, using different settings for the high-pass filter. The data is from GPS satellite 17, year 2014, DOY 50, site Honningsvag. (a), (b), and (c) have been calculated using a cutoff frequency of 0.2, 0.3, and 0.4 Hz, respectively.
(a)
(b)
(c)Taking into account the first and the second sensitivity analysis, the acceptable range for the computation ofσ
ϕ is f
c
∈
[
0.1,0.3
] and 〈
T
¯
〉
∈
[
30,60
]. These ranges give the maximum variations in the data.


## 7. Implementation and Analysis

### 7.1. Data Collection and Analysis

As mentioned in the introduction, the data used in this investigation are obtained from the Norwegian Ionospheric Scintillation Network. The reference stations used are encircled as shown in Figure1.For implementation and software packages used to compute the scintillation indicesS
4 and σ
ϕ, the interested reader is referred to Appendix C.


### 7.2. Interpretation of Results

The test results show that all detrending filtering techniques produce almost the same results. The main difference exists in statistical methods used in data editing. That is, how outliers are detected and removed, the missing observations, and glitch detection and removal.The new approach is robust against outliers, discontinuity, and missing observations. Figure14 shows how algorithms (FFT convolution, IIR, and local kernel regression) handle a discontinuity with a magnitude of 0.6 rad. The new approach handles perfectly the discontinuity and avoids generating a false scintillation due to discontinuity.Figure 14
The panels show 1 minute of phase data that contains a discontinuity, filtered using different techniques. (a) has been filtered using FFT convolution filter. (b) has been filtered using a IIR Butterworth filter. (c) has been filtered using the new approach.
(a)
(b)
(c)Figure 15
I95 index generated by the CPOS software for north of Norway.Figure 16
I95 index generated by the CPOS software for south of Norway.


### 7.3. Remarks

Some remarks and suggestions are as follows:(I)
In order to capture false scintillation due to low satellite elevation angles, the cutoff elevation angle was set to zero. Figure5 shows the events which are presented by shadowed rectangle.
(II)
In this paper, we have computed scintillation indices by means of different detrending and filtering techniques. Identifying the weaknesses and the strengths of each method, we can compute a reliable index that can be used in further analysis, for instance, to carry out the classification of the level of the ionosphere disturbances.
(III)
Nonparametric local regression with optimal smoothing parameters can be used as scintillation indices computation.
(IV)
As can be seen in Figure14, the new approach handles discontinuities very well, while the other filter techniques generate artificially high σ
ϕ values. The new approach computes σ
ϕ
=
0.0748 rad while the other algorithms compute values of 0.226 and 0.212 rad, respectively. Using the traditional techniques this discontinuity would be falsely reported as a scintillation event.


## 7.1. Data Collection and Analysis

As mentioned in the introduction, the data used in this investigation are obtained from the Norwegian Ionospheric Scintillation Network. The reference stations used are encircled as shown in Figure1.For implementation and software packages used to compute the scintillation indicesS
4 and σ
ϕ, the interested reader is referred to Appendix C.


## 7.2. Interpretation of Results

The test results show that all detrending filtering techniques produce almost the same results. The main difference exists in statistical methods used in data editing. That is, how outliers are detected and removed, the missing observations, and glitch detection and removal.The new approach is robust against outliers, discontinuity, and missing observations. Figure14 shows how algorithms (FFT convolution, IIR, and local kernel regression) handle a discontinuity with a magnitude of 0.6 rad. The new approach handles perfectly the discontinuity and avoids generating a false scintillation due to discontinuity.Figure 14
The panels show 1 minute of phase data that contains a discontinuity, filtered using different techniques. (a) has been filtered using FFT convolution filter. (b) has been filtered using a IIR Butterworth filter. (c) has been filtered using the new approach.
(a)
(b)
(c)Figure 15
I95 index generated by the CPOS software for north of Norway.Figure 16
I95 index generated by the CPOS software for south of Norway.


## 7.3. Remarks

Some remarks and suggestions are as follows:(I)
In order to capture false scintillation due to low satellite elevation angles, the cutoff elevation angle was set to zero. Figure5 shows the events which are presented by shadowed rectangle.
(II)
In this paper, we have computed scintillation indices by means of different detrending and filtering techniques. Identifying the weaknesses and the strengths of each method, we can compute a reliable index that can be used in further analysis, for instance, to carry out the classification of the level of the ionosphere disturbances.
(III)
Nonparametric local regression with optimal smoothing parameters can be used as scintillation indices computation.
(IV)
As can be seen in Figure14, the new approach handles discontinuities very well, while the other filter techniques generate artificially high σ
ϕ values. The new approach computes σ
ϕ
=
0.0748 rad while the other algorithms compute values of 0.226 and 0.212 rad, respectively. Using the traditional techniques this discontinuity would be falsely reported as a scintillation event.


## 8. Conclusion

We have shown that the new approach based on local regression with kernel smoothing and with AIC/AICC for bandwidth parameter selection can be used for computation of scintillation indices for high latitude data. The new approach can be analytically related to the existing filtering methods and is shown to produce highly correlated results with the traditional approaches. However, the new approach shows superior handling of discontinuities.We have shown that applying the right detrending and filtering techniques to the scintillation data at high latitude one can obtain reliable scintillation indices (S
4 and σ
ϕ). Clearly, the studies show that the phase scintillations are dominant for these data sets.For the derivation ofS
4, the digital filter FIR (Butterworth, Chebyshev, and Elliptic) and FFT convolution used to implement the low-pass filter work well with minor differences. For σ
ϕ, the digital filter IIR (Butterworth and Elliptic) and FFT convolution work well. The difference exists in statistical methods used to compute standard deviations and how outliers, glitches, and missing observations are handled. Poor methods can bias the estimation process. In addition, the study shows that the derivation of scintillation indices S
4 and σ
ϕ is sensitive to averaging interval and the cutoff frequency. Figure 5 shows the false scintillation for S
4. In this case, all elevation angles are used under computation.For the data sets classified as high ionospheric activity, we have defined the empirical distribution of the scintillation indicesS
4 and σ
ϕ for the Norwegian Regional Ionospheric Scintillation Network, located at high latitude (61.9
9
∘N to 70.9
8
∘N). Heavy tailed Frechet/log-normal distribution (B.5) is a good model for phase scintillation σ
ϕ index as Figure 8 confirms. The distribution of S
4 follows the Nakagami distribution, given by (B.1) and shown by Figure 6.



---

*Source: 3582176-2016-08-07.xml*