Challenges of Real-Time Monitoring of Ionospheric Perturbations and TEC Fluctuations with GPS Single Station

1. Introduction

Ionosphere has a big influence on degradation of accuracy and reliability of the positioning with trans-ionospheric radio signals. Its influence is very critical. In the post processing mode, ionospheric characteristics can be determined very easily, but in near- and real-time it is a very challenging task. Determination of ionospheric characteristics requires permanent monitoring in real time, and nowcasting and forecasting of ionospheric indices. For now, we are familiar with the well known ionospheric indices:

• number of electrons along the signal propagation path: Total Electron Content ($TEC$);

• rapid and fast fluctuation of Radio Frequency (RF) signals’ amplitude and phase ($S4$ and ${\sigma }_{\varphi }$) and

• rate of change of $TEC$ ($ROT$).

In order to obtain information about the state of the ionosphere using single station GPS observations, we are developing and constantly upgrading our iono-tools module that is a part of the in-house academic software TUB-NavSolutions.

Previously, we presented some of the possibilities and methods to monitor ionospheric amplitude scintillation [1] and now we are dealing with $TEC$ calculations, smoothing and levelling methods. As our algorithms and software are being tested, $TEC$ estimation performance has been analysed in a simulated real-time mode. The achieved results are described here.

$TEC$ values can be calculated from code-or carrier phase measurements. Usage of the carrier phases requires challenging ambiguity fixing, while $TEC$ derived from code-phases are noisy.

Thus, for monitoring of the TEC in real-time the decision has been met to smooth code $TEC$ using the carrier $TEC$. An approach, the so called levelling of $TEC$ derived from the carrier phases using the code-phases, has been applied to overcome challenging fixing of carrier-phase ambiguity terms to integer numbers. Selected methods of smoothing of code $TEC$s (levelling of carrier-phases) are described and compared with TEC available online at the Center for Orbit Determination in Europe (CODE) in the IONosphere Map Exchange Format (IONEX). In this paper some results of the analysis of applicability and performance of the algorithms, using GPS observations from two selected days have been presented:

• 15^th July 2012, day of the year 197 (DOY 2012:197) and

• 16^th July 2012, day of the year 198 (DOY 2012:198).

Data that has been taken into analysis has been collected at the Kiruna station in the polar area.

To track intensity of Earth’s geomagnetic activity and to detect geomagnetic storms, two indices have been chosen: Ap and Kp. It is worth to mention that even though geomagnetic and ionospheric storms are related, geomagnetic storms refer to disturbance of Earth's magnetic field, and ionospheric storm is a disturbance of the ionosphere [8].

2. Methodology

2.1. Total Electron Content estimation from GNSS single station measurements

Slant $TEC$ has been estimated directly from GNSS dual frequency carrier- and code-phase in-situ measurements using the following equations:

$${TEC}_{carrier}=\frac{1}{40.28}\frac{f_1^2{f}_2^2}{f_2^2-f_1^2}\left({\varphi }_2-{\varphi }_1\right)$$E1

$${TEC}_{code}=\mathrm{ }\frac{1}{40.28}\frac{f_1^2{f}_2^2}{f_2^2-f_1^2}\left(P_1-P_2-c\frac{f_2^2-f_1^2}{f_2^2}{T}_{gd}\right)$$E2

Where

• ${\varphi }_1$ and ${\varphi }_2$are carrier-phase observations $\left[m\right]$,

• $P_1$ and $P_2$ are code-phase observations $\left[m\right]$,

• $f_1$ and $f_2$ are frequencies of $L1$ and $L2$ signals $\left[Hz\right]$,

• $T_{gd}$ is transmitter’s Estimated Group Delay between $P_1$ and $P_2$ measurements $\left[s\right]$,

• $c$ is velocity of the $[m/s]$.

$TEC$ derived from the carrier-phases are ambiguous and those ones derived from code-phases are very noisy. To take advantage that the ambiguous carrier-phases measurements have low noise and that the noisy code-phase measurements are not ambiguous, some smoothing and levelling methods have been applied.

2.2. TEC smoothing and levelling methods

$TEC$ is calculated from both carrier- and code-phases using equations 1 and 2. After the both are calculated, levelling offset has been added according to the following formula:

$${TEC}_{slant}={TEC}_{carrier}+{TEC}_{offset}$$E3

Where value of the offset depends on the method has been applied.

2.2.1. Method I

Calculation of the levelling offset between code- and carrier-phases derived $TEC$ along the whole arc according to [7].

$${TEC}_{offset}=\frac{1}{N}\sum _{i=1}^N{TEC}_{code}\left(t_i\right)-\mathrm{ }{TEC}_{carrier}\left(t_i\right)$$E4

Where $t_i$ indicates time period in which $TEC$ has been observed and $N$ is number of calculated $TEC$ in the whole arc of satellite visibility.

2.2.2. Method II

Levelling approach suggested by Jakowski [5]. Here the offset is defined by the following formulas:

$${TEC}_{offset}={TEC}_{code}-{TEC}_{carrier}$$E5

$${TEC}_{offset,\mathrm{ }i+1}=\mathrm{ }\frac{i}{1+i}{⟨{TEC}_{offset}⟩}_i+\mathrm{ }\frac{1}{i+1}{TEC}_{offset,\mathrm{ }i+1}$$E6

2.2.3. Method III

Smoothing and filtering of code- with carrier-phases according to Hatch algorithm [3]. $TEC$ derived from code-phases are being smoothed by previous $TEC$ derived from carrier-phases observations.

$$\widehat{{TEC}_{code}\left(t_i\right)=\frac{i}{N}+\frac{N-1}{N}}\left[\widehat{{TEC}_{code}}\left(t_i-1\right)+\mathrm{ }\left({TEC}_{carrier,i}-{TEC}_{code,i}\right)\right]$$E7

For real-time usage, when method III will be applied (Hatch filter), length of the filter must be reinitialized whenever cycle slip is detected. Detection of the cycle slips must be performed in all smoothing approaches.

2.3. Amplitude scintillation

For calculation of the amplitude scintillation, the $S4$ index was derived from signal power calculated from the $I$ (in-phase) and $Q$ (quadrature) components of bandpass signal (eq. 8):

$$s\left(t\right)=I\left(t\right)+jQ\mathrm{ }\left(t\right)$$E8

Power of the signal is now derived as:

$$P\left(t\right)=\frac{I^2+Q^2}{2}$$E9

The scintillation index has been calculated using following formula:

$$S4=\sqrt{\frac{⟨P^2\mathrm{ }\left(N\right)⟩-{⟨P\mathrm{ }\left(N\right)⟩}^2}{{⟨P\mathrm{ }\left(N\right)⟩}^2}}$$E10

where

$$⟨P\mathrm{ }\left(N\right)⟩=\frac{1}{N}\sum _{i=0}^{i=N}P\mathrm{ }\left(i\right)$$E11

and

$$⟨P^2\mathrm{ }\left(N\right)⟩=\frac{1}{N}\sum _{i=0}^{i=N}{P}^2\mathrm{ }\left(i\right)$$E12

The argument $i$ indicates the observation epoch number within the data time span, and $N$ is the number of available observations.

According to the formula given by A.J. Van Dierendonck [2] $S4$ index has been calibrated by the influence of ambient (thermal) noise.

$${S4N}_0=\sqrt{\frac{100}{SNR}\left[1+\frac{500}{19SNR}\right]}$$E13

Now, the final calculation of $S4$ index is done by:

$${S4}_{calibrated}=\sqrt{{S4}^2-{S4N}_0^2}$$E14

2.4. Rate of change of Total Electron Content

In order to trace ionospheric irregularities [11] and to provide spatial variation of electron density [9], the rate of change of $TEC$ ($ROT$) is introduced. The equation below (eq.13) describes estimation of the $ROT$ parameter:

$$ROT=\frac{{TEC}^{t+i}-{TEC}^t}{\left(t+i\right)-t}=\frac{\left({\varphi }_1^{t+i}-{\varphi }_2^{t+1}\right)\left(\mathrm{ }\frac{1}{40.28}\mathrm{ }\frac{f_1^2{f}_2^2}{f_1^2-f_2^2}\right)-\mathrm{ }({\varphi }_1^t-{\varphi }_2^t)\mathrm{ }\left(\mathrm{ }\frac{1}{40.28}\frac{f_1^2{f}_2^2}{f_1^2-f_2^2}\right)}{\left(t+i\right)-t}$$E15

In our case we are calculating between epochs changes of $TEC$ every one second. That gives us simplified equation of $ROT$ without denominator because it is always equal to $1 sec$ (assuming there are no gaps).

3. Case study

Data processed for this analysis were collected at the Kiruna station in Sweden ($67.5026°N$, $20.2437°E$) with approximate position depicted in the Figure 1. In Kiruna there is a GNSS continuously operating station of the German Aerospace Center (DLR), division in Neustrelitz. The station is working for the Space Weather Application Center-Ionosphere (SWACI). It is configured for ionospheric scintillation monitoring and the observables (code- and carrier phases, and the $I$- and $Q$-amplitudes) are recorded with $50 Hz$ sampling rate.

Figure 1.

Position of Kiruna station in Sweden

Raw values of the observables collected at the station have been processed with university software being under development at Technische Universität Berlin (TUB) TUB-NavSolutions. The software module iono-tools is depicted in the Figure 2.

Figure 2.

Basic architecture of TUB-NavSolutions software iono-tools module

Data taken into the analysis is from days with high geomagnetic activities. In the Figure 3. $Ap$ and $Kp$ indices are shown for only couple of days: 14-17. July 2012. In this paper we will focus on 15^th and 16^th July.

$TEC$ values derived directly from the observations were compared with slant $TEC$ derived from Global Ionospheric Maps (GIM) of the Center for Orbit Determination in Europe (CODE), Astronomical Institute of the University of Bern following Schaer [10].

Figure3.

Daily averaged indices of geomagnetic activity

4. Results

Figure 4 and 5 show $TEC$ for both days and for two selected satellites: PRN 09 and PRN 27. The light blue curve represents slant $TEC$ values derived from CODE GIM and the magenta curve represents elevation angle. In the both Figures slant $TEC$ derived from GIM ($TEC$ from IONEX) is represented by a smooth curve. However, in the both Figures sudden and small peaks, in the all three smoothed time series of $TEC$ are visible. These anomalies tell us that most probably some ionospheric disturbances appeared at that time. $TEC$ derived from observations, smoothed and levelled with all three methods coincide with each other but all of them deviate from the GIM ones for more than $5 \left[TECU\right]$.

In the Figure 4 interesting are two deviations both detected roughly between 9:00^h and 11:00^h UTC (marked with green circles). This is approximately the same time when $Ap$ and $Kp$ indices reached values that indicate a strong geomagnetic storm.

Figure 4.

$TEC$ calculated and interpolated for 15^th and 16^th July 2012 for satellite PRN 09

Figure 5.

$TEC$ calculated and interpolated for 15^th and 16^th July 2012 for satellite PRN 27

$ROT$ values calculated from observations indicate disturbed ionospheric conditions for both days. In the both Figures, 6 and 7, we can notice sudden peaks and variations which confirm that at that time ionospheric perturbations took place.

If we compare $TEC$ from the left panel of the Figure 4 with $ROT$ from the left panel of the Figure 6 it is seen that larger oscillations appear in both time series at the same time. Special warning (based on strong $ROT$ variations) comes a little bit before 11:00^h.

Figure 6.

$ROT$ for 15^th and 16^th July for satellite number PRN 09

Figure 7.

$ROT$ for 15^th and 16^th July for satellite number PRN 27

Similar behaviour and coherence is visible also in other time series of $TEC$ and $ROT$ (marked with red circles in Figures 4 and 5)

• in the right panels of the Figures 4 and 6 a little bit after 9:00^h and a little bit after 10:00^h,

• in the left panels of the Figures 5 and 7 with big peaks a little bit before 9:00^h, around 9:30^h and before 11:00^h,

• in the right panels of the Figures 5 and 7 with sudden and big oscillation a little bit after 9:00^h and constant oscillations between 10:00^h and 11:00^h.

In the Figure 8 ionospheric amplitude scintillation parameter ($S4$) has been depicted. $S4$ values are displayed for each day and for all satellites in view. Even though ionosphere amplitude scintillation is less intense in polar regions, a few higher values (above $0.6$) may be observed indicating ionospheric perturbations.

Figure 8.

Amplitude scintillation index $S4$for 15^th and 16^th July for all satellites

5. Discussions and conclusions

The data from only two selected days has been post processed. Both days, 15^th and 16^th July 2012, have been selected on the base of geomagnetic indices. There are clearly visible similarities between time series of $TEC$ and $ROT$ values derived from observations, on one side, and geomagnetic indices ($Ap$ and $Kp$), on the other side. Variations in the time series of $ROT$ are very similar to those in the time series of $EC$. Even small jumps (peaks), visible in $TEC$ plots, can be assumed as correlated with $ROT$ oscillations.

Some $TEC$ smoothing and levelling methods have been tested here in order to select the most appropriate one for our real -and near-real time applications. All three tested methods give very similar results of the final $TEC$ values. It has been found that the method $II$ and $III$ fulfil requirements for usage in real-time. In the Fig. 4 and 5 results from the two methods, $II$ (green line) and $III$ (blue line) are displayed. Both curves are overlapping. That is why we can have impression that blue curve does not exist, but looking for numerical values it is seen that they differ between each other on second place after decimal point only.

There are easy seen biases between $TEC$ derived from our obervations and interpolated using CODE. A source of the biases is not identified yet because not enough data was available up to now. Investigation of it will be continued.

The above described draft results of investigations allow to assume that $TEC$ and $ROT$ variability can be used for detection of perturbed ionospheric conditions and probably for issue of warnings for real-time users.

Investigation of applicability of $TEC$and $ROT$for real-time warnings on ionospheric perturbations will be continued using GPS data collected at stations located in equatorial, mid-latitude and polar areas.

Acknowledgments

The research is funded by the FP7 People Programme through the Marie Curie Initial Training Network TRANSMIT-Training Research and Applications Network to Support the Mitigation of Ionospheric Threats.

We would like to thank Space Weather Application Center Ionosphere (SWACI) operated by German Aerospace Center (DLR) in Neustrelitz for the SWACI data sets.