1. Introduction
The modernisation of the Global Navigation Satellite System (GNSS) and the constant improvement of their performance is driving the design of a plethora of applications with strict requirements on accuracy, precision and reliability of the estimated position [
1]. Despite the impressive improvement of the quality of the broadcast signals and new processing techniques, the propagation of the signal through the ionosphere is still one of the limiting factors for precise positioning. In fact, it is well known that the GNSS signals while travelling through the upper layers of the athmosphere are affected in terms of delay, amplitude and phase, which turn out to be different with respect to the assumed propagation in free-space. Compensation of the ionospheric delay is applied by now in any GNSS receiver, either using models or external assistance for single frequency receivers, and it can be compensated by means of the iono-free combination of the measured pseudoranges in multi-frequency receivers [
2]. However, due to the presence of a non-uniform distribution of the electron density in some layers of the ionosphere, electromagnetic signals are scattered in random directions, thus reaching the Earth through multiple paths. Signals on each path add in a phase-wise sense and interfere with each other, so that the received signal is then affected by deep amplitude fading and phase fluctuations.
GNSS signals are indeed affected by scintillation effects, and their quality at the receiver might be strongly degraded. While moderate scintillations only increase noise in code and phase measurements, strong scintillation events severely disrupt the capabilities of GNSS receivers. Amplitude fading and rapid phase changes impact the receivers’ tracking loops capabilities to follow the signal dynamics, resulting in cycle slips and Loss-of-Lock (LoL) in the tracking stage of the receiver [
3]. Phase range measurements errors can be from 3 to 10 times larger [
4], whereas, in the most severe situations, when more than one satellite-receiver link is concerned and the satellite visibility is poor, the receiver might completely fail [
5]. At the same time, thanks to their global coverage and continuous availability, GNSS signals provide an excellent means to monitor the activity of the ionosphere. With multiple satellites from multiple constellations, the physical effects can be simultaneously measured through many pierce points. GNSS receivers have then become a convenient, cheap and widespread tool to study such phenomena [
6,
7].
In this paper, we propose a technique to assess the intensity of scintillation on the signal amplitude, which is an alternative to the classical methods that require a GNSS receiver employing closed-loop architectures. Rather then relying on the traditional commercial receiver architecture, in which only the correlator outputs are available to the end user, this approach fully exploits the flexibility of a receiver architecture entirely designed according to software-radio approach [
8]. Samples of the received signal at Radio Frequency (RF) are collected, sampled and digitized, at sampling frequencies of the order of millions of samples per second. Then, the bit-stream is stored on memories and processed by a fully-software GNSS receiver. This approach has gained consensus in recent years also in the space weather community [
9,
10,
11], and it allows for an easier implementation of novel algorithms, specifically tailored to the scintillation monitoring and mitigation [
12]. As a consequence, such techniques can be robust by design to the strong scintillations events causing LoLs in classical receivers’ architectures.
The paper is organized as follows:
Section 2 will briefly recall the well known parameters used for scintillation detection and classification.
Section 3 will discuss the flexibility given by a software-defined-radio approach for the design of novel methods to process the digitized GNSS signal, including the possibility to undertake different approaches for the scintillation monitoring.
Section 4 discusses the impact of scintillation on the classical closed-loop architectures. The core of the paper is introduced in
Section 5, where an open loop architecture is proposed to implement a robust interference monitoring platform, and in
Section 6 new metrics are presented. Finally,
Section 7 will present the results obtained processing real data collected during scintillation events.
2. GNSS Scintillation Measurement
As previously recalled, ionospheric scintillations are rapid and hard-to-predict fluctuations of the amplitude and phase of radio signals propagating through the atmosphere. The propagation through ionosphere is one of the main error sources affecting the GNSS signal processing stages of the receiver. In fact, the signal is potentially affected by strong degradation, introducing significant problems in any GNSS application requiring high accuracy, precision, reliability and continuity. On the other hand, scintillation mitigation is still an open issue in most of the commercial devices, mainly due to the inherently random nature of such events, which makes modelling impossible [
13]. Scintillations are then, at the moment, one of the main limiting factors for high accuracy applications.
Scintillation events are traditionally characterized by means of two indices evaluated over an 60 s observation interval.
The amplitude scintillation index, named , estimates the amount of amplitude fluctuations. Such an index is computed as the normalized standard deviation of the detrended Signal Intensity (SI), derived from in-phase (I) and quadrature (Q) prompt correlation outputs, normalized by its mean value.
The phase scintillation index denoted measuring the impact on the phase of the signal; it is equivalent to the standard deviation of the detrended carrier phase measure.
The majority of the works related to scintillation detection and monitoring are based on the analysis of the value of these two indices, often comparing them with predefined thresholds. However, the estimation of
and
requires the GNSS receiver to be in lock conditions, i.e., the values of I and Q correlators must be available. In the case of strong scintillation events, receivers designed according to the classical closed-loop architecture might not be robust enough to keep the correct tracking of the signals, incurring in a LoL event. If this happens, the receiver is not able to provide the values needed to characterize the the scintillation event, which is, in a generic way, classified as “strong”. At the same time, the computation of
and
requires a detrending operation, which envisages the use non-trivial mathematical operations. Although this process has been standardized and is largely used, several authors pointed out that it can introduce artifacts in the signal processing chain falsifying the scintillation estimation [
14,
15]. This opens the way to the definition of alternative measures of scintillation [
14,
16,
17,
18].
Some of these methods start from the observation that, in the presence of scintillations, the statistics of the samples of the amplitudes change over time and it is correlated to the values taken by the
index. In statistical terms, according to [
4], scintillation-derived amplitude fluctuations follow a Nakagami-m distribution, the shape of which depends on the severity of the fluctuations. This is depicted in
Figure 1, where the histograms of the samples for different
values are displayed. The method introduced in this paper aims at estimating the statistic of the amplitude samples in the presence of scintillations, building on top of that a new metric for the characterisation of the scintillation event, as will be explained in the following sections.
4. Impact of Strong Scintillations on GNSS Receivers
One of the consequences of scintillation is LoL of the tracking stages of the receiver, affecting the quality of navigation solution and also making impossible the calculation of the scintillation activity through the classical indices. In case of LoL, the computation of indices using tracking loop outputs is not reliable, thus leading to false alarms and missed detections of the monitoring unit.
An example of LoL due to strong amplitude scintillation in closed loops is reported in
Figure 3. From top to bottom, the values of the I and Q prompt correlators, the estimate of the
, the values of the Phase Lock Loop (PLL) discriminator output, the estimate of the Doppler frequency
and the estimated
value are reported.
Figure 3a reports the tracking results obtained running a classical second order tracking loop. After about five minutes, the fading induced by amplitude scintillation is so strong that the PLL loses the lock on the signal. As a consequence, the
values provided by the scintillation monitoring engine, and computed upon the values of the I and Q correlators, cannot be trusted. As a benchmark, the results obtained running an experimental third order loop, especially configured in terms of bandwidth and loop parameters to be able to track such extreme fading situations, are reported in
Figure 3b. The third order loop is able to keep the lock on the signal, even though the quality of the measurements is low. In this case, the estimate of the
can then be considered valid. It is worth noting that, for the first five minutes, the two architectures provide the same value of the
index, while, after the loss of lock, the values provided by the 2nd order PLL are false, and, in particular, they underestimate the intensity of the scintillation event.
In general, under significant scintillation events, cycle slips, degradation of
and corruption of data bits affects the quality of the estimated position (since some satellites might have to be excluded from the solution). In addition, the estimations of
values for the affected channels are wrong, as shown in
Figure 3a.
Various techniques are proposed in the literature to reduce the probability of a LoL in the presence of scintillation. The straightforward solution is indeed to exploit a third order PLL [
26,
27], as shown by the previous example, but also advanced receiver architectures have been proposed, [
12]. However, the only solution able to completely overcome the problem of LoLs is to abandon the closed loop implementation in favour of an open loop architecture [
18,
28]. It is, of course, well known that closed loop solutions are optimal in terms of limiting the jitter on the code and phase estimations, but they show limitations for the specific use in scintillation monitoring receivers in the presence of strong events.
In the following sections, an open loop architecture will be exploited to extract from the incoming signal the information related to the amplitude scintillation intensity, defining a metric alternative to the
parameter. This architecture has been firstly presented in a previous paper of the authors [
18], but has been refined and extended in the present work. A similar approach has been proposed in other works (such as [
29,
30]); however, the architecture proposed in this study overcomes several limitations, among which the need of a fine time assistance.
5. Open Loop Processing Architecture
The key concept of GNSS open loop processing, which is depicted in
Figure 4, is that, as long as the receiver is static, its reference position and time are known, and a network connection providing A-GNSS messages is available; some of the components of the signal can be considered as known, or properly predicted. The knowledge of such components can then be exploited to remove them from the incoming signal, leaving only the random contributions that are due to the thermal noise and other disturbances, among which are the scintillation effects. Scintillation contributions can then be separated from the GNSS signal, which can be considered a sort of deterministic, periodic contribution in the overall received signal.
Let
be the discrete version of the received signal at the output of the analog-to-digital converter of the front-end which operates at a sampling frequency
on the GNSS signal at IF.
is the number of satellites in view, each of one contributing to the overall received signal; the pedix
i refers to the
i-th satellite;
is the navigation message;
is the code delay;
is the spreading code;
is the carrier frequency, which includes the nominal intermediate frequency
and the Doppler frequency shift
;
is a generic phase contribution; and
is thermal noise.
The architecture does not include either a PLL or a Frequency Lock Loop (FLL) demodulating the signal to baseband by removing the residual carrier at frequency
. The residual carrier is wiped-off using the Doppler information from assistance. Such assistance information is a good estimate of
, but it does not include the spurious effect that can be due to the instability of the front-end clock. For such a reason, a good quality clock with negligible frequency drift has to be used to drive both the sampling stage and the frequency references in the front-end. The effect of the non-perfect carrier removal and the issues related to the Doppler shift estimation will be discussed in
Section 5.2. Let us assume, for sake of simplicity at this stage, that the estimation of the Doppler frequency is perfect and then that the carrier the wipe-off demodulates the signal to baseband. The mixing generates images of the spectrum at a harmonics frequency, which are removed by the accumulation process, acts as a low pass filter. The signal
, after carrier wipe-off, and assuming perfect time synchronization and ideal filtering, can then be expressed as:
and is, in general, a complex number, due to the multiplication times a complex exponential.
As a following step, if a data channel is considered, the navigation message has to be removed, (e.g., exploiting assistance message and the reference time at the receiver). Small errors in this process can be tolerated, as they introduce negligible losses, as it will be shown later. The signal
, after navigation message wipe-off, and assuming perfect time synchronization, can then be expressed as:
Because of the presence of thermal noise
, the Pseudo Random Noise (PRN) spreading code
is not visible and therefore it is not yet possible to estimate and remove it. In principle, the PRN code could be removed by a wipe-off operation, as in a pure open loop receiver, but this would require sub chip-level time synchronization and fine assistance [
31,
32]. This condition is difficult to be satisfied only exploiting proper A-GNSS information and good knowledge of the position of the receiver. Therefore, a different strategy for code removal has to be identified, based on a blind self-estimation of the code. Since thermal noise can be modeled as a zero mean Gaussian random variable, an accumulation process can increase the SNR (Signal-to-Noise Ratio) of the signal up to the point in which the code chips emerge from the noise floor. The digital signal after carrier and data wipe-off,
, is divided in chunks of length
, where
is the spreading code periodicity, as, for example, 1 ms for the GPS C/A signal. Supposing a period is composed of
N samples,
L periods are coherently summed together, obtaining the signal of length
N
as depicted in
Figure 5.
The scope of this coherent accumulation is the estimation of the GNSS spreading code, and this technique has been widely used in the past for the blind identification of the codes transmitted by new GNSS satellites [
33] or to detect anomalies and evil waveforms [
34,
35]. The estimated binary GNSS code can be extracted by means of the sign operator:
The spreading code is then wiped-off from the GNSS signal, by multiplying the binary sequence
times the signal at output of the accumulation,
, thus obtaining a complex signal
affected only by a residual noise component and by scintillation:
The digital signal is eventually fed to the scintillation metric block.
The code retrieval process is depicted in
Figure 6: few chips of the the signal after accumulation,
, are depicted in blue, while the estimated binary spreading code
is depicted in red. On the reconstructed code
, some residual distortion can be noted due to any other un-modeled effect that cannot be averaged to zero by the accumulation process, and, as well, due to some residual noise contribution.
5.1. Trade-Off between Noise Average and Observation Rate
The signal will be processed by the scintillation analysis block that will perform a statistical analysis of the samples. For this reason, the noise contribution still affecting the signal after the accumulation has to be sufficiently small, so that it does not mask or bias the results.
Given a fixed sampling frequency , increasing L, the noise contribution tends to be averaged out, but this would require many segments of received signal, i.e., a large observation window . Due to the non-stationarity of the scintillation process, the statistical distribution of the signal amplitude changes over time. This limits the increasing of L that, if chosen to be too large, would not catch the significant changes on the statistic distribution since they would be all included in the same observation window. The parameter L must trade off this resolution in time of the analysis with the capability of reducing the noise effect.
As a trade-off between the noise averaging and the resolution in time for the GPS C/A code, , corresponding to five seconds of the received signal, was heuristically chosen. Therefore, from now on, signals obtained accumulating over 5000 periods of ms of signal, are adopted for scintillation analysis of the GPS C/A code.
5.2. Transition Samples
In order to be able to coherently accumulate the code periods, the carrier removal must be very effective. The Doppler frequency shift affecting the signal can indeed be derived from A-GNSS information, combining the information on the satellites’ position and velocities and the known user position. Therefore, the residual carrier
can be ideally wiped-off. Similarly, the navigation message can be perfectly removed, by exploiting A-GNSS aiding. However, the Doppler effect does not only change the carrier frequency, but also the code rate, causing it to stretch and compress. Supposing
to be the nominal chip rate,
the transmission frequency, it can be found that the actual chip rate at time
t,
is:
where
is Doppler frequency shift on the carrier at time
t.
Given a fixed sampling frequency , the time varying code rate induces a change in the number of samples per code chip from period to period. Furthermore, and are not likely to be multiples of each other. The result is that the same chip, in subsequent PRN code periods, might not be represented by the same number of samples. In other words, subsequent sampled chips might not be perfectly aligned, thus threatening the effects of the coherent accumulation of the code periods.
Ideally, after the accumulation process, an exact representation of the spreading code would be obtained, as depicted by perfect rectangular wave plotted in red in
Figure 6. However, in a realistic situation, as shown in blue in
Figure 6, the non-perfect alignment of the samples leads to “transition samples” at the boundaries of the chips that are wrongly reconstructed by the quantization performed by the sign function.
Transition samples cannot be avoided, unless a very low value of
L is used, which could, in turn, prevent a correct extraction of the code from the noise. The effect of transition samples over the overall statistic of samples is depicted in
Figure 7, where histograms of
with and without transition samples are reported.
Figure 7a clearly shows a tail in the distribution towards zero, due to the transition samples.
If transition samples are not removed from , a residual modulation remains in after the code wipe-off. Their presence can indeed alter the statistic analysis performed to extract scintillation metrics. For this reason, a simple algorithm to discard transition samples has been implemented, discarding both the sample before and after a sign transition. This strategy avoids the definition of thresholds of a multilevel non-uniform quantizer for the removal of the transition samples, which could introduce further artifacts in the signal processing chain.
6. Scintillation Statistical Analysis
The scintillation analysis is performed through a goodness of fit process taking into account the statistical distribution of the samples in the cases with and without scintillation effects.
A scintillation-free reference curve is created best-fitting a Gaussian curve to the distribution of the samples in absence of scintillation. corresponds to the distribution of the samples for the time interval under test. If fit is accurate, is a good approximation, then there is no scintillation, and, vice versa, scintillation is declared present.
The model curve for estimating
is
in which
a,
b,
c are estimated on the basis of an iterative method that minimizes the sum
S of squared residuals
where
Figure 8 reports an example, depicting the reference curve
and the histogram of a set of samples that includes both scintallated and non-scintillated samples. As it can be noted, the Gaussian fitting is able to match only the non-scintillated part of the histogram, while in the presence of scintillation, the samples are distributed according to the typical Nakagami-m probability density function, thus deviating from
obtained in case of no scintillation.
In statistics, to evaluate the goodness of fit, the parameter
is used, where
is the mean of
, and
are the values of the best-fitting Gaussian distribution for the non-scintillated case.
can vary in the interval
. The greater
, the better the fitting and vice versa. A measure of diversity can then be defined through a novel parameter denoted
and defined as
7. Validation and Use of the Open-Loop Metrics
In this section, in order to validate the definition of the the new metrics and demonstrate their use, real datasets collected in low-latitude regions are used. A first dataset is derived from the data collection campaign collected in the frame of the MImOSA project [
36] by means of a GNSS SDR data acquisition system [
8]. Data are related to GPS L1 C/A signal collected along equatorial regions at −22.120045° N, −51.408671° E on 25 March 2015. The signal was sampled at 5 Msamples/s, with a resolution of 16 bits, and demodulated at baseband. To measure the effects of scintillations, the
and
parameters are employed and compared. The
has been computed by means of a software routine processing the I and Q correlator outputs provided by the closed loop tracking of a GNSS software receiver, and following the methodology explained in
Section 3.
has been derived following computations of
Section 6. In such a data collection, several satellite-receiver links are affected by scintillations.
Figure 9 plots the
and
obtained for PRN1 and PRN3. In particular,
values are shown in blue and with a rate of
Hz, whereas, in red,
R4 values are reported with a rate of
Hz.
It can be noted that
and
are somehow correlated. However,
is obtained at a higher rate than
, and, being based on the open-loop architecture, it does not require the tracking to be locked, i.e., it can be obtained also in case of strong scintillation events that threaten the correct operation of a GNSS receiver, as detailed in
Section 7.1.
In order to demonstrate that
is actually carrying the same information of
, a filtering operation is applied to
values to resemble
values. A moving average window of width
and with amplitude
is used. Parameters
and
were obtained minimizing the root mean square error between the two curves over a large dataset of scintillation events representative of all possible dynamics. Similarity between the curves in
Figure 10 shows how
follows the same trend of the classical
index and thus it carries the same kind of information. This comparison validates the definition of
as an alternative to
that can be obtained at a higher rate. When, in case of strong events, the estimation of
is not available or not reliable,
can still be calculated and it extrapolates the same kind of information carried by
also in these threatening scenarios.
Furthermore, many recent works on scintillation are questioning the validity of the classical algorithm for
estimation, mainly because of the detrending procedures involved, which, if not properly tuned, might mask or falsify scintillation events. A different index, independent from
, not requiring detrending operations, and exploiting technological advances such as software defined receivers, might help in this process of redefining the metrics for the characterization of the scintillation events [
14,
37].
7.1. Case of Loss of Lock
In order to test the goodness of the open loop approach in the presence of a loss of lock, the same case study reported in
Section 4 is analyzed. It has to be remarked how
, being based on an open loop architecture, is not prone to LoL, and it can work also in case of very strong scintillations.
Figure 11 reports the value of the
filtered metric, obtained by processing the raw signal samples by means of the open loop receiver. This is compared to the values of the
obtained by the second order closed loop receiver, as reported in
Figure 3a, where a LoL occurred after about five minutes. In order to benchmark the results, a reference
, obtained by a commercial scintillation monitoring receiver, is plotted, corresponding to the
estimate obtained running a robust third order loop, as shown in
Figure 3b.
As already observed, the estimate from the closed loop cannot be trusted after minute 5, due to the incorrect values of the I and Q correlator outputs. On the contrary, the metric derived from the open loop processing follows the same trend of the reference , thus confirming the validity of the proposed technique in the presence of strong amplitude scintillations inducing losses of lock in traditional receiver architectures.