EDP Sciences logo
Submit your paper
Search
Menu
All issues Volume 14 (2024) J. Space Weather Space Clim., 14 (2024) 17 Full HTML
Issue
J. Space Weather Space Clim.
Volume 14, 2024
Topical Issue - Observing, modelling and forecasting TIDs and mitigating their impact on technology
Article Number 17
Number of page(s) 15
DOI https://doi.org/10.1051/swsc/2024017
Published online 01 August 2024

© M. Guerra et al., Published by EDP Sciences 2024

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Travelling Ionospheric Disturbances (TIDs) are plasma density fluctuations propagating as waves through the ionosphere at a wide range of frequencies and velocities (Hooke, 1968; Hines, 1974; Kersley & Hughes, 1989; Hocke & Schlegel, 1996; Hernández-Pajares et al., 2012; Habarulema et al., 2022) and are usually classified according to those parameters into Small, Medium, and Large scale TIDs (SSTID, MSTID, and LSTID, respectively). Specifically, LSTIDs are characterized by periods between 45 and 180 min and phase velocities that are in the 300–1000 m/s range (Hocke & Schlegel, 1996; Borries et al., 2009; Katamzi et al., 2012; Thaganyana et al., 2022). MSTID instead fall in the 10–60 min period and 50–350 m/s phase velocity range, respectively (Ogawa et al., 1987; Hernández-Pajares et al., 2012). Finally, SSTIDs have not been studied as extensively as LSTIDs and MSTIDs, but their wavelengths are below 100 km (Tsybulya & Jakowski, 2005; Yin et al., 2019; Boyde et al., 2022). TIDs are mainly generated by the coupling between the neutral and ionized atmosphere (Hines, 1974; Hocke & Schlegel, 1996; Zakharenkova et al., 2016; Astafyeva, 2019), with gravity, acoustic, and acoustic-gravity neutral waves forcing a movement on ions due to collisions with neutral particles. Therefore, all phenomena generating atmospheric neutral waves can be detected by their ionospheric signature, provided that the wave amplitude is large enough and that they propagate up to such heights. Generally, LSTIDs are related to increased geomagnetic activity, which heats the neutral atmosphere due to Joule heating, Lorentz forces, and particle precipitation within the auroral oval and thus generates equatorward propagating gravity waves (Hunsucker, 1982; Jing & Hunsucker, 1993; Borries et al., 2009; Zakharenkova et al., 2016; Cesaroni et al., 2017; Frissell et al., 2022). MSTIDs, instead, are related to natural hazards, such as earthquakes (Leonard & Barnes Jr., 1965; Astafyeva et al., 2009; Astafyeva & Heki, 2009), tsunamis (Artru et al., 2005; Zhang & Tang, 2015), volcanic eruptions (Themens et al., 2022; Verhulst et al., 2022, Madonia et al., 2023), tropospheric storms (Lay et al., 2013; Chou et al., 2017), and other phenomena caused by changes in the ionizing radiation from the Sun, e.g., solar terminator passages (Song et al., 2013), solar eclipses (Chimonas & Hines, 1970; Zhang et al., 2021), solar flares (Zhang et al., 2019), and electrodynamical processes, such as Perkins instability associated with E- and F-region coupling processes (see, e.g., Makela & Otsuka, 2012). Therefore, due to the strict relationship between natural hazards and TIDs, it could be possible to leverage ionospheric measurements to issue alerts related to natural hazards (Ravanelli et al., 2021) and to increase our understanding of the processes behind the Lithosphere-Atmosphere-Ionosphere coupling (LAIC).

In addition to being of scientific interest, TIDs are also a driver for GNSS and HF degradation. For the mid-latitude areas, where other phenomena (e.g., plasma irregularities) are less frequent than at low and high latitudes (Basu et al., 2002; Kotulak et al., 2020), TIDs are the main cause of HF and GNSS disruption. LSTID can be responsible for a variation in TEC up to 2 TEC units (1 TECu = 1016 electrons) for mid-latitude regions (Borries et al., 2009) and up to 3–4 TECu for low-latitude ones (Valladares et al., 2009; Habarulema et al., 2016). MSTIDs are usually responsible for variations of tenths of TECu, but they can reach up to a few TECu during solar maximum (Hernández-Pajares et al., 2006, 2017) and they can also cause an increase in cycle slip occurrence (Poniatowski & Nykiel, 2020). TIDs are sometimes regarded as a “silent killer of accuracy” or “short-range catastrophe” (Fung et al., 2022) as they are hardly detectable when they affect the GNSS signals received on the ground. TIDs are considered a major nuisance for any system using trans-ionospheric radio wave propagation in the mid-latitude regions (Belehaki et al., 2020). Specifically, Timoté et al. (2020) show that MSTID activity causes an increase in the positioning error for Network Real-Time Kinematic (NRTK), while Poniatowski & Nykiel (2020) show the same impact on precise point positioning (PPP). It is thus essential to detect and characterize the whole plethora of TIDs in real-time, making it possible to develop mitigation techniques (such as the one presented by Hernández-Pajares et al., 2017), forecasting algorithms and natural hazard early warning techniques. A review of the widely adopted methods for TID detection, with a focus on real-time operations, can be found in Belehaki et al. (2020).

Many techniques are available to extract periodic features from time series and have been extensively used in previous studies. These techniques are usually divided into two categories: filtering and detrending techniques. In this paper, we considered four detrending and two filtering techniques. The first detrending technique taken into consideration is the largely used moving average/running mean, which has been extensively exploited for the detection of both MSTIDs and LSTIDs due to its extreme simplicity and computational lightness (Nishioka et al., 2013; Cherniak & Zakharenkova, 2018; Ferreira et al., 2020; Otsuka et al., 2013, 2021). The second detrending technique is the double time difference, which subtracts from each value an average value of the previous and posterior measurements. DD has been used mainly for MSTID studies due to its uncompared simplicity and its ability to amplify the features falling in the period band of interest (Hernández-Pajares et al., 2012, 2017; Zhang & Tang, 2015). The last two detrending techniques are based on polynomial fitting: polynomial detrending and Sovitsky-Golay filtering (Schafer, 2011). Polynomial detrending works by removing the best-fit polynomial from the time series; in the literature, polynomials of different orders have been used to detect both MSTIDs and LSTIDs (Astafyeva & Heki, 2009; Pradipta et al., 2014, 2016; Habarulema et al., 2016, 2022). The Savitsky-Golay filter is a more sophisticated version of the polynomial detrending where the polynomial fitting is performed on a sliding window (Wang et al., 2007; Katamzi et al., 2012; Coster et al., 2017). Katamzi et al. (2012) used the Savitsky-Golay filter because “it is adaptive and gives compartively good perfomance” to the moving average and is more robust against boundary effects and non-linearities. Moving to filtering techniques, the two algorithms considered are fast iterative filtering and finite impulse response band-pass filtering. FIF will be briefly explained in the data and methods section since it is a novel technique used in a few geophysical studies related to TIDs (Verhulst et al., 2022). FIR-based BUTF instead has been leveraged in many studies to highlight periodic features falling in the frequency band of interest (Calais et al., 2003; Wang et al., 2007; Komjathy et al., 2012; Van De Kamp et al., 2014; Bukowski et al., 2024). Usually, FIR-based bandpass is employed because of its ability to highlight features in a specific frequency band while stopping those falling outside of the pass band.

Moving to RT and NRT, there are currently only a few GNSS-based algorithms for the detection of TIDs related to natural hazards, such as VARION (Savastano et al., 2017; Savastano & Ravanelli, 2019; Ravanelli et al., 2021) and GUARDIAN (Martire, Krishnamoorthy, Vergados, Romans, Szilágyi, Meng, Anderson, Komjáthy and Bar-Sever, 2022). The VARION algorithm was originally developed to detect tsunamis in NRT but is now used to detect several kinds of natural hazards (Ravanelli et al., 2023). It is based on the single time difference of the carrier phase measurements geometry-free linear combination (GFLC): this represents the VARION core observation (dTEC/dt). Thanks to the time differentiation, all the constant biases (ambiguity and inter-frequency biases) are removed and can be used in real time. The GFLC time difference can be further filtered or detrended using the most common techniques: Butterworth band-pass (Ravanelli et al., 2023) and 8th-order polynomial filter (Savastano et al., 2017) were widely explored. The TID detection module of the GUARDIAN algorithm instead works by applying a 4th-order Butterworth high-pass filter to the uncalibrated (i.e., including the biases) GFLC of phase measurements, giving near-real-time TEC estimates.

Given the importance of TID detection and characterization, it is necessary to know the reliability of the techniques used to extract the wavy feature of interest from the TEC time series. Maletckii et al. (2020) have already performed a similar study, focusing on the errors (amplitude and cross-correlation) induced by the different techniques. They investigate the performance of different detrending and filtering techniques. In contrast, this study does not stop at the impact of filtering techniques: it investigates how specific parameters (height of the IPP and elevation cut-off) and different NRT observables affect the extracted signal period and amplitude. To our knowledge, no previous study tried to quantify the impact of the Doppler effect on the extracted period given by the relative movement of the TID and the IPP. Finally, it is important to stress that all currently available NRT TID detection algorithms tend to focus mainly on the wave period characterization while disregarding the precise estimation of the wave amplitude, which might turn very useful in mitigation techniques and for natural hazard studies (Astafyeva, 2019; Astafyeva et al., 2022).

The paper is structured as follows: Section 2 introduces the data used to generate our database of synthetic arcs. Then, it presents the different detrending techniques and explains how we quantified the impact of such techniques, ionospheric parameters and NRT observable on the extracted wave parameters. Section 3 covers the statistical analysis results. Finally, a brief discussion on the possibility of generalization of the findings of this work to other geographical areas is performed, together with some discussion on the reliability of the assumptions.

2 Data and methods

The following section covers the preparation of synthetic TEC arcs and the analysis performed on the obtained database of GNSS measurements. We decided to leverage calibrated TEC arcs (Ciraolo et al., 2007; Cesaroni et al., 2015, 2021) obtained from a GNSS station in Southern Italy to obtain data representative of the real ionospheric background ionization and IPP movement. Once the database of calibrated verticalized TEC is obtained, some pre-preprocessing (smoothing to remove high-frequency components) is performed before adding the modelled TIDs to better compare the different detrending techniques. After this, two wave fields are applied separately to the pre-processed vertical TEC. Now, having a database with many different observation arcs, it is possible to apply all the detrending techniques and compare their output with the known ground truth (the wave field reconstructed at the IPP).

2.1 Arc reconstruction

The choice of a synthetic database is driven by the necessity of a known ground truth (our reconstructed wave field shown in equation (3). That is why having a smoothed real ionospheric background with a superimposed modelled wavy feature has been chosen as a compromise to obtain our two synthetic ionospheric datasets (one for the MSTID and one for the LSTID case). Since we decided to use real measurement as the background, we considered two years of data (2009, solar min, and 2014, solar max). The two periods defined by different solar activity are chosen to ensure the presence of different background ionization and day-to-day variability. As the amplitude of our waves is fixed, having different values of background ionization will better represent the real ionosphere. It will thus test detrending techniques in a more realistic scenario.

All the raw GNSS data were gathered by a geodetic station belonging to the INGV Rete Integrata Nazionale GNSS (RING, ring.gm.ingv.it; RING WORKING GROUP, 2016), which is located close to Altamura (AMUR: 40.91°N, 16.60°E, see Figure 1), Apulia, Southern Italy. The RING network has been proven to provide high-quality RINEX data to monitor TEC over Italy (Cesaroni et al., 2021), and a real-time mapping over the region is currently available in the electronic Space Weather upper atmosphere (eswua.ingv.it) data portal managed by INGV (Cesaroni et al., 2020). The data used were in the form of 30s RINEX files, and calibrated TEC was retrieved through the Ciraolo calibration technique (Ciraolo et al., 2007; Cesaroni et al., 2021; Tornatore et al., 2021). Such calibration technique leverages the simultaneous observations at different frequencies, which allows obtaining an observable free of errors that are not a function of the carrier frequency. This is performed, arc by arc, by levelling the carrier phase GFLC to the GFLC of code measurements. The observable obtained is the following:

L̃ARC=sTEC+bR+bS+〈εP〉,$$ {\mathop{L}\limits^\tilde}_{{ARC}}={sTEC}+{b}_R+{b}_S+\left\langle {\epsilon }_P\rangle, $$(1)

thumbnail Figure 1

Reconstructed LSTID wave field calculated over Europe on a 0.5×0.5 lat/lon grid. The left panel represents the wave field calculated at t equal to 0, while the right panel represents the wave field at t equal to a half period. The black dot represents the location of the RING station AMUR.

where 𝑏𝑅 and 𝑏S are the receiver and satellite Differential Code Biases (DCBs), and sTEC is the TEC along the slant path connecting the given receiver and satellite. As one can notice from equation (1), the new observable is not ambiguous, and as suggested by Braasch M. (1996) we can disregard carrier phase noise and multipath. Following Ciraolo et al. (2007) we merge 〈𝜀𝑃〉 with DCBs, obtaining a new error, 𝛽𝐴𝑅𝐶, as shown in equation (2):

L̃ARC=sTEC+βARC,$$ {\mathop{L}\limits^\tilde}_{{ARC}}={sTEC}+{\beta }_{{ARC}} $$(2)

A thin shell ionosphere at a defined height is assumed to proceed with the Ciraolo calibration process. In this way, it is possible to derive a mapping function that relates the sTEC to the vTEC and depends on the height of the ionospheric thin shell layer and satellite elevation (Mannucci et al., 1998). Now, thanks to the thin-shell ionosphere assumption, it is possible to derive 𝛽𝐴𝑅𝐶 for each receiver-satellite pair. Once the calibration and mapping process is completed, a 15000 real vTEC arcs database is available. The arc database includes only GPS satellites, and it is available at https://zenodo.org/doi/10.5281/zenodo.10074542.

As both filtering and detrending techniques are compared, the high-frequency components are removed from the arc database by performing a Gaussian-weighted moving average. Filtering techniques can remove the high-frequency components, while the detrending techniques cannot: this pre-processing should ensure that the lower frequencies induce no bias in the amplitude error.

For our study, two TID wave fields are considered; thus, we will have two different arc databases related to LSTID and MSTID. The two wave field parameters are set according to climatological values for the European sector for both MSTIDs (Hernández-Pajares et al., 2006, 2012) and LSTIDs (Hocke & Schlegel, 1996; Borries et al., 2009; Thaganyana et al., 2022), as shown in Table 1.

Table 1

Wave field parameters (amplitude, wavelength, azimuth, speed, and period) for the two different TID scenarios.

Figure 1 shows a graphical representation of the LSTID wave model calculated on a 0.5° × 0.5° latitude-longitude grid. The blue and yellow bands represent the plasma depletion/increase caused by the wave field. The two panels show phase values equal to 0 and π. As the considered TIDs are defined by different periods, the Gaussian-weighted moving window size used for pre-processing is set according to 1.33 times the given TID period (22.5 and 100 min for the MSTID and LSTID scenarios, respectively). A 2D plane wave field, presented in equation (3), is used to add the TID to the smoothed background arc databases. We assume that the wave field propagates at a given height equal to the IPP height (HIPP). This assumption is justified by the ionospheric plasma conditions at mid-latitudes, where most of the ionized matter is localized in a thin layer; thus, we assume that the TID propagates in the thin ionospheric layer (Mannucci et al., 1998). Moreover, HIPP is usually set by the user according to realistic values of the height of the F2 ionospheric layer (Astafyeva et al., 2022), which further justifies our assumption. We then model the TEC variation (dTEC) over space (x, y) and time (t) induced by a wave like the following:

dTEC(x,y,t)=Asin(2πxλcosθ+2πyλsinθ+2πvtλ ,)$$ {dTEC}\left(x,y,t\right)={Asin}\left(\frac{2{\pi x}}{\lambda }{cos\theta }+\frac{2{\pi y}}{\lambda }{sin\theta }+\frac{2{\pi vt}}{\lambda }\enspace \right) $$(3)

where A is the wave amplitude, λ is its wavelength, θ is the azimuth of propagation, v is the TID speed, x and y are latitude and longitude converted to km of arc over the Earth’s surface.

2.2 Detrending techniques

Once the arc database is generated, the following step is to define the detrending techniques that are applied to all the wavy TEC arcs obtained. To extract wavy features from the vTEC signal, we applied the six different detrending/filtering techniques introduced in the first section to the reconstructed arcs:

Moving average.The multi-order numerical difference.Third-order Savitzky-Golay filter.8th-order polynomial detrending.Finite Impulse response band-pass filter.Fast Iterative Filtering band-pass filter.

While techniques 1–5 are well-known and widely applied in many scientific fields, the FIF-based band-pass filter requires further explanation. FIF is a novel signal decomposition technique capable of decomposing non-stationary and non-linear signals into simple oscillatory (around zero) components called intrinsic mode functions (IMFs) (Cicone & Zhou, 2021), as shown in equation (4):

s(t)=∑i=1NIMFIMF(t,ν)i+res.$$ s(t)=\sum_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{IMF}}}\mathrm{IMF}{\left(\mathrm{t},\mathrm{\nu }\right)}_{\mathrm{i}}+\mathrm{res} $$(4)

In equation (4), (𝑡) is the signal to decompose, 𝑁𝐼𝑀F is the total number of IMFs, ν is the given IMF frequency at time t, and res corresponds to non-oscillatory residuals, which is assumed to be the background trend. To leverage FIF as a band-pass filter, we need to sum all IMFs defined by frequencies inside the range of interest at every time step. FIF is a technique similar to Empirical Mode Decomposition (Huang et al., 1998), with the advantage of being around a hundred times faster. FIF has proven to be particularly suitable for analyzing ionospheric data, to reveal multi-scale features of the ionospheric medium (see, e.g., Ghobadi et al., 2020; Spogli et al., 2021, 2019; Urbar et al., 2022) and signatures of TID due to LAIC events (Verhulst et al., 2022). As FIF is based on the Fast Fourier Transform, it assumes periodicity at the boundaries of the signal under study. Stallone et al. (2020) proposed an extension approach that drastically reduces boundary errors without introducing spurious modes in the decomposition. Thus, the signal extension increases the signal length with a half-sample antisymmetric extension to avoid the boundary effect induced by the decomposition technique.

The main parameters of the six techniques were set according to a preliminary study that focused on extracting the best parameters for those techniques that do not rely only on the period band or sliding window size (SGOLAY and POLY). The period bands were instead set according to the usual LSTIDs and MSTIDs characteristics, with the most impactful being the period. The period band of BUTF and FIF for MSTIDs was set between 10 and 40 min, while the 45–90 min band was used for LSTIDs. The polynomial degree was set to 10th for the MSTID scenario and 5th for the LSTID one. The MA was performed within a sliding window of 60 min for the LSTID scenario and 30 min for the MSTID one. The DD time separation for the MSTID scenario was set according to the value (300s) used in GNSS derived TEC based climatological studies on MSTIDs that exploited the DD detrending technique (Hernández-Pajares et al., 2006, 2012). For the LSTID scenario, as no example of DD application for the large-scale case was found, we set the time separation (30 min) equal to half the size of the MA window. Finally, from the preliminary study, we found that the best SGOLAY performance for the LSTID and MSTID scenarios was obtained by fitting a 2nd order polynomial in a sliding window twice the size of the MA one (30 and 60 min).

Figure 2 shows the reconstruction and detrending process for a single arc under the MSTID propagation scenario. Specifically, panels a, b and c of Figure 2 show an example of the different steps of the reconstruction process for the MSTID scenario. Calibrated vTEC, obtained from the Ciraolo calibration (Figure 2a), is the real ionospheric background. Instead, Figure 2b and 2c are the smoothed background vTEC before applying the wave field, and the sum between the latter and the TID model is calculated at each IPP by applying equation (3). Figure 2d shows the detrended (FIF-based band-pass) version of the synthetic arc along with the TID model calculated at the IPP locations.

thumbnail Figure 2

The top panels (a, b, c) show the different arcs of TEC considered in the reconstruction process under the MSTID scenario. Specifically, panel (a) shows the output of the GG calibration algorithm (vTEC), panel (b) shows the moving average smoothed vTEC arc and panel (c) shows the smoothed arc summed with the TID wave model (blue line of panel (d)). The bottom panel (d) instead shows the detrended (FIF-based band-pass) version of the synthetic arc (orange line) along with the TID wave model calculated at the IPP locations (blue line).

2.3 Detrending techniques comparison

Once the two synthetic arcs databases had been detrended with the six techniques mentioned in the previous section, it was possible to investigate the amplitude error (AME) and time domain error (TDE) induced by the wave extraction technique. Note that TDE refers to time-shifting or contraction/extension in time of the decomposed signal, an artefact that prevents accurately deriving the TID speed from interferometry-like techniques. For the time domain error (TDE), we used one minus the normalized cross-correlation at lag 0. According to this, AME and TDE are expressed as:

AME=xr-x$$ {AME}={x}_r-x $$(5a)

TDE=1-xr∙x (xr∙xr)(x∙x),$$ {TDE}=1-\frac{{x}_r\bullet x\enspace }{\sqrt{{(x}_r\bullet {x}_r)(x\bullet x)}} $$(5b)

where x represents the detrended arc and 𝑥r the reconstructed wave field. Because of the different nature of the two errors, while the number of AMEs is equal to the sample number, there is only one TDE for each arc. Finally, we looked at the error distribution and statistics of AME and TDE to better investigate the general picture of the error.

2.4 Doppler effect and induced period error

Once AME and TDE caused by the different detrending algorithms have been ranked, the focus is moved to the effect on the wave period induced by the underlying nature of the measurement. As said, we assumed that the TID propagates in a thin shell layer; thus, different heights of such a shell will be related to different IPP speeds (Hernández-Pajares et al., 2006). Since the IPPs are moving across the TID, the measured wave will suffer a Doppler effect directly related to the IPP speed and, thus, to the satellite elevation and height of the ionospheric shell (e.g., IPP velocity increases for lower elevation angles and higher shell heights). Savastano et al. (2017) show that, by setting an incorrect ionospheric height, the TID estimated speed differs from the real one. Despite that, TID speed is usually derived from interferometry techniques (Afraimovich et al., 1998; Maletckii and Astafyeva, 2021), and thus, the main issue would be the difficulty in correlating waves that show different periods due to the Doppler effect. Detrending techniques are unnecessary to investigate this phenomenon because the Doppler effect depends only on the relative velocity between the IPP and the TID. Therefore, since the TID velocity vector is set and remembering that we assumed the TID to propagate precisely at the height of the ionospheric shell, we only need to derive the IPP velocity vector for each data sample. Once the IPP velocity vector is evaluated starting from the consecutive IPP coordinates, we can calculate the Doppler effect thanks to the following equation:

f'=f(1+vREL|v⃗TID|) where: vREL=v⃗TID∙v⃗IPP|v⃗TID.|$$ {f}^{\prime}=f\left(1+\frac{{v}_{{REL}}}{\left|{\vec{v}}_{{TID}}\right|}\right)\enspace {where}:\enspace {v}_{{REL}}=\frac{{\vec{v}}_{{TID}}\bullet {\vec{v}}_{{IPP}}}{\left|{\vec{v}}_{{TID}}\right|} $$(6)

In equation (6), vREL is the TID-IPP relative velocity, vTID is the TID velocity vector, f is the original wave frequency and 𝑓′ is the apparent frequency. As the IPP velocity strongly depends on the satellite elevation and IPP height, we investigated the Doppler-induced error statistics for different values of the two parameters. Finally, it is possible to infer the expected error given the ionospheric shell height and the elevation cut-off angle through statistics.

2.5 Near real-time derivable observables and related amplitude error

In this section, we present two different NRT-derivable observables, and we propose an evolution of one of the two that should ensure a higher level of accuracy in the estimation of the extracted signal amplitude. As of now, the first two observables have been used extensively for NRT and post-processing applications, and many examples are available in the literature (Astafyeva et al., 2022; Brissaud & Astafyeva, 2022; Ravanelli et al., 2021; Maletckii & Astafyeva, 2021; Martire et al., 2023). To evaluate the performance of such observables, a statistical analysis is performed by comparing the signal extracted from the three NRT observables with the one derived from a widely used post-processing technique proposed by Ciraolo et al. (2007) The first NRT-derivable observable we have considered is the sTEC derived from phase measurements, which can be retrieved in real-time through the Networked Transport of RTCM via Internet Protocol (NTRIP). To compute the sTEC starting from phase measurements, we need to apply the following equation:

sTECPH=1A*f12*f22f12-f22*(L1*λ1-L2*λ_2),$$ {sTE}{C}_{{PH}}=\frac{1}{A}*\frac{{f}_1^2*{f}_2^2}{{f}_1^2-{f}_2^2}*({L}_1*{\lambda }_1-{L}_2{*\lambda }\_2) $$(7)

where 𝐴 = 40.308 m3 s−2, L1 and L2 are the phase measurements, and 𝜆1 and 𝜆2 are the wavelength of the two frequency bands 𝑓1 and 𝑓2. Once the phase-derived sTEC is obtained, the TEC value of the first data point is subtracted from the entire time series to limit the unknown bias characteristic of phase measurements. By detrending this phase-derived sTEC, it is possible to extract the wavy features correctly in the time domain. However, its reliability in extracting amplitudes is not comparable due to the slanted nature of the observable. Therefore, one must verticalize it to extract the amplitude reliably (especially for low elevations where there is the least accordance in the sTEC and vTEC magnitude). The second observable is thus obtained by converting sTECPH to vTECPH through the mapping function:

vTECPH=sTECPH*cos⁡(arcsin⁡(RecosθRe+HIPP)),$$ {vTE}{C}_{{PH}}={sTE}{C}_{{PH}}*\mathrm{cos}(\mathrm{arcsin}(\frac{{R}_e{cos\theta }}{{R}_e+{H}_{{IPP}}})) $$(8)

Where Re is the Earth’s radius, θ is the elevation angle, and HIPP is the ionospheric shell height. Note here that when applying the mapping function to sTECPH, the bias will not be constant anymore due to the dependence on the elevation angle. By looking at the mapping function, it is possible to see that the induced amplitude error will be minimal for high elevation. That is why we propose the third and novel approach, where an initial guess is performed by estimating the initial value of the sTECPH time series through an ionospheric model. One of the main constraints of the ionospheric model is that it must be extremely quick (hundredths of seconds per TEC arc) in computing the initial bias. This constraint is extremely important for NRT applications. Once the initial guess is gathered from the model, the sTEC is verticalized by applying equation (8); in principle, thanks to this approach, the amount of error induced by the verticalization of the initial bias should be lower, with a minimum for high elevations and a maximum for low ones (this observable will be called vTECNe). The ionospheric model we used is NeQuick2 (Nava et al., 2008), which was selected based on its computational lightness and speed at calculating the TEC (some ms per arc), given the receiver and satellite coordinates. Figure 3 shows four different observables belonging to the same arc, three of which have been obtained from the three different NRT techniques under the MSTID scenario, and the fourth is the calibrated vTEC, which, due to its post-processing nature, is considered as the ground truth. Focusing on the NRT observables, the blue dashed line represents sTECPH. Following, the blue dotted line is vTECPH (so, the red one times the mapping function), while the dash-dotted blue line is vTECNe. At first glance, one can notice that the detrended arcs (Figure 3b) are almost identical for the values close to the middle arc time. This divergence is due to the elevation (orange line in Figure 3a) decreasing while reaching the arc boundaries since every arc starts when the satellite rises and ends when the satellite disappears below the horizon. Moreover, it is worth noticing that Figure 3b shows that the wave period changes throughout the arc duration, with the first half decreasing while the second half increases. Both aspects will be covered in a more generalized way in the results section.

thumbnail Figure 3

The top panel (a) shows the different NRT-like observables studied together with the satellite elevation. The solid blue line is an example of a reconstructed arc, just like the one in Figure 2c. The orange line instead is the satellite elevation. The bottom panel (b) focuses instead on the detrended version of the four observables presented in panel (a), which share the same colour code. We zoomed into dTEC peaks corresponding to low, medium, and high elevations to better visualize the slight differences between the detrended observables.

As for the Doppler error section, we are not investigating the reliability of the six detrending techniques but rather which NRT observables best resemble the results obtainable in a post-processing study. Therefore, the detrending technique is fixed for both scenarios. The technique used is SGOLAY because it proved to be the most accurate in the detrending techniques comparison analysis. Once the NRT arcs are detrended, the AME is computed. Note here that the different NRT approaches do not cause different effects on the TID period. However, they amplify or attenuate the signal amplitude depending on the sign of the mapped bias. That is why the NRT error analysis only covers the AME while disregarding the TDE.

3 Results

3.1 Detrending techniques comparison

First, the amplitude error induced by the different detrending techniques is investigated. Figure 4 shows the overall AME statistics for the six techniques. From the Cumulative Distribution Functions (CDFs) in the top panels, one can quantify how likely it is to expect an AME error smaller than the one of interest. For example, it is possible to know that 80% of the time, the AME is smaller than the corresponding x value for the technique of interest by fixing a cumulative P(x) of 80%. Looking at the CDF under the LSTID scenario and fixing a cumulative probability of 80%, one can observe that FIF and SGOLAY are the best techniques AME-wise, with an x value of 0.125 TECu (around 35% of the original amplitude). In contrast, the numerical difference technique (DD) performs the worst (almost 0.5 TECu, more than 100% of the LSTID amplitude). This pattern is confirmed by the bottom-left panel, where one can see that FIF and SGOLAY show the smallest percentiles. Now, focusing on the MSTID scenario (top-right panel), we can note that the general picture is similar to the LSTID one, with FIF and SGOLAY being again the most reliable techniques (80% of the AME are around 0.05 TECu, 25% of the original MSTID amplitude). Instead, DD is the least accurate (0.1 TECu, 50% of initial amplitude). Regarding the other three techniques, one can remark that none clearly outperforms the others. In particular, the blue and purple lines for the LSTID scenario of Figure 4 almost perfectly overlap. Before looking into the correlation error, we would like to stress that the DD technique, which is the one performing worst, is not meant to extract the amplitude reliably since it is supposed to amplify up to 2× the wavy feature defined by the period set as the main technique parameter (time separation).

thumbnail Figure 4

Amplitude error statistics for the six techniques considered. MSTID (right column) and LSTID (left column) parameters correspond to the ones presented in Table 1; the ionospheric thin shell height was set to 350 km. The top panels show the CDFs of AME for the different techniques (each line corresponds to the CDF of a given technique). X-axis shows the AME error in TECu, while the Y-axis corresponds to the probability P(AME) of having an error equal to or smaller than the corresponding X value. The bottom panels, instead, present the median along with their related percentiles (5, 16, 84, and 95) for the different filtering/detrending techniques.

The violin plots shown in Figure 5 are used to investigate the correlation error statistics. The violin plots effectively present the probability distribution and the overall statistics (mean, median, and interquartile) in a compact way. First, focus on the LSTID scenario (left panel). As for the AME, FIF and SGOLAY perform best in preserving the original wavy feature in the frequency domain, with FIF showing slightly better performance. The mean and median for the FIF distribution (cyan) are the smallest, and the same applies to the interquartiles. In addition, thanks to the violin plot, it is possible to observe that FIF has the highest amount of correlation errors close to zero. In addition, FIF also shows the thinnest distribution for high correlation error values, meaning it is the least prone to large frequency distortions. Before investigating the MSTID case, we want to stress that every distribution is highly asymmetrical, confirmed by the significant difference between mean and median. Focusing on numerical values, FIF, MA, and DD show a comparable correlation error mean under the MSTID scenario. In contrast, in the LSTID case, FIF and SGOLAY outperform the other techniques.

thumbnail Figure 5

Violin plot for the correlation coefficient error. The X-axis shows the correlation error, while each violin (colour-coded according to different techniques) is a box plot where the shape represents the probability distribution of the correlation error for the given technique. The black diamond and line are, respectively, the median and mean of the correlation error distribution, while the shaded area represents the interquartile range. The left panel represents the LSTID scenario, while the right one is for MSTID.

Considering both TID scenarios, we can conclude that FIF and SGOLAY are the best techniques again, ranking first and second for both the LSTID and MSTID scenarios.

To conclude, considering the findings derived from Figures 4 and 5, we can conclude that, among the techniques considered, FIF and SGOLAY have proven to be the best at extracting wavy features from post-processing observables, assuring the best performance in estimating the amplitude while reliably preserving the wave frequency. In addition, concerning the computational speed of the techniques, SGOLAY would be the technique of choice since it is around 300 times faster than FIF. Nevertheless, some considerations about the DD technique must be made. As previously stated, the DD technique is designed to amplify features up to 2×, with the highest degree of amplification for features defined by a frequency equal to the one corresponding to the moving window. Thus, including the DD technique in the AME ranking is somehow unfair because it is specifically designed to act on the amplitude of the feature of interest.

3.2 Doppler effect and induced period error

After ranking the reliability of the six detrending techniques in estimating the wavy features (amplitude and period), the focus is moved to the Doppler-induced period error results. Figure 6 shows that the CDFs of the period error change depending on the TID scenario, IPP height, and elevation cut-off. The CDFs of period error for the LSTID scenario and an IPP height of 250 km are shown on the top left. The first clear result is that the lower the elevation cut-off, the higher the expected period error on average. Specifically, the 80% cumulative P(×) corresponds to a period error of 9, 11, and 15% for 60-, 40- and 20-degree elevation cut-off, respectively. Now, we can compare the scenario (LSTID, 250 km) with the bottom-left panel, which is the same plot but for an increased IPP height (350 km). Thanks to this comparison, one can quantify how much the increased IPP height influences the extracted period.

thumbnail Figure 6

Cumulative distribution functions (CDF) of the Doppler error for different TID scenarios, heights of the IPP, and elevation masks. Plots on the right/left correspond to MSTID/LSTID, while top/bottom plots are for an IPP height of 250/350 km. The differently coloured lines show the CDF for different elevation masks, where blue, red, yellow, purple, and green are 20, 30, 40, 50, and 60 degrees, respectively. The X-axis represents the percentage period error, while the Y-axis shows the cumulative probability.

Specifically, considering the 80% cumulative P(×), we obtain period errors of 20, 14, and 12% of the original period for 20, 40, and 60 degrees, respectively. As expected from theory and confirmed by the statistics, a higher IPP height causes the period error to increase on average.

Now, focus on the MSTID scenario. The first clear observation is that the period error almost doubles on average compared to the LSTID case. This behaviour was strongly expected since the Doppler-induced period error is directly connected to the TID speed, as seen in equation (6); therefore, the slower the TID, the higher the Doppler effect. Numerically speaking, the same 80% cumulative P(×) related to the same elevation cut-offs as for LSTID are 21, 22.5, and 33% for an IPP height of 250 km and 32, 34, and 53% for 350 km under the propagating MSTID scenario. Thus, the slower MSTID shows period errors that are doubled with respect to the faster TID, a ratio that corresponds to the two TIDs speed ratio, 2 (See Table 1). Moreover, looking at the MSTID scenario, we can observe that the CDFs show a flex point, a feature that is not visible in the LSTID ones. A possible explanation for this flex point will be given in the Discussion section. Figure 6 allows us to quantify the expected period error for different realistic TID scenarios and geometrical parameters, such as IPP height and elevation cut-off.

3.3 NRT observable AME statistics

Lastly, the focus is moved to the amplitude error statistics when applying a fixed detrending technique to NRT available observables. Figure 7 shows the median/percentiles of the three observable CDFs for the two TID cases. Looking at the CDFs, it is possible to remark that vTECNe outperforms the other two observables in both the MSTID and LSTID scenarios. Focusing only on the LSTID CDFs, one can notice that sTECPH is way less accurate than the others. Specifically, the 80% cumulative P(×) amplitude errors for vTECNe is almost negligible (0.015 TECu, comparable to the noise level for phase measurements), while for STECPH and vTECPH it is respectively around 0.16 and 0.065.

thumbnail Figure 7

Amplitude error statistics for the NRT observables derived dTEC. The top panels show the CDF of the amplitude error, with each line representing a different observable (blue, red, and yellow for sTECPH, vTECPH, vTECNe, respectively). The X-axis shows the amplitude error, while the Y-axis is the cumulative probability. The bottom panels show the median and the 18–84 and 5–95 percentiles for the different techniques. Here, the Y-axis represents the amplitude error in TECu. The plots on the left and right correspond to LSTID and MSTID, respectively.

This difference between the three observables is confirmed further by the median and percentiles shown in the bottom-left panel of Figure 6. Contrarily, the performance of the two verticalized observables is similar under the MSTID scenario. Numerically speaking, fixing at 80% the cumulative probability, we get an AME of 0.075, 0.005, and 0.0025 for sTECPH, vTECPH, vTECNe. Furthermore, by looking at percentiles, it is possible to notice that the 16–84 percentiles for vTECPH and vTECNe are almost identical (both around 0.001 TECu, symmetric distribution), while the ones for sTECPH is clearly higher (0.05 TECu). This difference means that sTECPH is way more prone to large error values, while vTECNe and vTECPH manage to reduce its impact/number. To conclude, considering both scenarios, we can confidently affirm that vTECNe is a reliable choice for extracting the TID amplitude in an NRT application. In addition, by comparing the vTECNe amplitude error performance with the one related to the best technique, FIF, it is possible to see that the error induced by working in an NRT scenario is smaller than the AME induced by the detrending process. Specifically, the vTECNe AME is 10% and 5% of the FIF error for the LSTID and MSTID scenarios, respectively.

4 Discussion

In this paper, we analyzed some of the most common errors impacting TID detection techniques for post-processing and NRT applications. We decided to use reconstructed arcs to obtain quantitatively and qualitatively reliable results. Nevertheless, one may argue that removing the high-frequency components might affect the numerical results. In principle, such simplification should not strongly affect the AME results because the high-frequency components are assumed to be random; thus, it should add to a null statistical impact considering the size of the database. Before considering the AME and TDE results, we must keep in mind that the techniques investigated are not all of the same nature. FIF, for example, works with a frequency band, while the easier MA needs only a window size. This single parameter means that the attenuation for the frequencies that are not of interest is less strong than for filtering techniques. That is why we decided to consider a monochromatic sine wave, which should limit the impact of the differences between filtering and detrending techniques.

Now, thanks to the analysis accomplished, we can affirm that among the considered techniques, FIF and SGOLAY are the most reliable overall in estimating the wave amplitude and period. By looking at the literature on detrending technique performances applied to TID studies (Maletckii et al., 2020; Osei-Poku et al., 2021), we can see that, despite the different techniques involved, our results are quite in accordance with other studies. We must remember that the GNSS station is located in southern Italy, which falls in the mid-latitude region, where the ionospheric background is quite different from other latitudinal sectors (Li et al., 2022; Mansoori et al., 2013; Venkatesh et al., 2011). Thus, the results shown might not be numerically reliable for high- and low-latitude. Despite that, the general rank gathered from the AME, and TDE analysis should still be valid.

Following, we tried to quantify the impact of the observational geometry in the form of shell height and elevation mask. If we consider TIDs generated by tsunamis and earthquakes, the height of maximum coupling is lower than the F peak, and thus it is usually below 250 km (Astafyeva and Shults, 2019; Thomas et al., 2018). This should ensure minimal period error when such natural hazards generate TID. Moreover, the effect of the IPP movement does not affect the arrival time of the TID, but it only affects its duration. In contrast, the elevation mask is normally set according to the user’s needs and habits and is normally higher than 20°. The exception is the GNSS-based tsunami detection (Martire et al., 2022; Savastano et al., 2017), where the lower the elevation, the further away from the coast we are sensing. However, lowering such a cut-off below 20° might induce multipath and other artefacts. Now, thanks to the Doppler error results, we can know in advance what degree of period error to expect, given our experimental setup (station latitude). Nevertheless, it should be noted that we have only studied two TID cases; thus, the exact numerical values will change if the TID azimuth and velocity differ from the considered ones. To conclude, the strange shape of the MSTID CDFs visible in Figure 6 is likely due to the nature of the IPP trajectories, which are directly related to the GNSS orbital parameters. Due to the MSTID speed being comparable to the IPP one and because of the azimuth being fixed, there might be a significant difference in the error statistics for different orbital planes. Thus, if we consider the CDF as a sum of the error CDF for each orbital plane, one might be largely shifted towards larger errors, causing such flex points. The same does not apply to LSTID because the wave speed is always larger than the IPP one, limiting the degree of error induced. Thanks to this analysis, we now know that the Doppler error can be easily as dramatic as 50% of the original period for MSTID detected at low elevations and high HIPP. In addition, we must point out that for TID velocity retrieval, there are techniques that should mitigate the Doppler effect, which spoils the quality of measurements (Penney and Jackson-Booth, 2015).

The final section showed that the proposed NRT observables performed, and some interesting aspects deserve further discussion. As presented in the results section, the proposed vTECNe observable showed an extremely low amplitude error, especially under the MSTID scenario. Specifically, 95% of the errors are lower than the usual noise level for 30s differential phase measurements, around 0.05 TECu (Coster et al., 2012). Thus, we can confidently affirm that leveraging such an observable in an NRT scenario would provide performance comparable to the post-processing ones. Moreover, the time needed to obtain the vTECNe observable is comparable with the ones related to the other two observables; thus, one can obtain a better degree of precision in exchange for a negligible increase in computational time. Since we only studied 30s data, the same might not be true for 1s measurements, as the typical noise level on 1s data is lower than in 30s ones. As for the Doppler effect, the numerical values might not be accurate for equatorial areas because the ionospheric background is less homogeneous than at mid-latitude (Eastes et al., 2019). In addition, due to the equatorial ionization anomaly, the reliability of NeQuick at low latitudes is not as good as at mid and high latitudes (Angrisano et al., 2013). Despite that, even if the error is likely bigger, we expect it not to be higher than the noise level, and the general rank should not change.

5 Conclusions

The main conclusions of this study are the following:

  • Among the techniques considered, FIF and SGOLAY proved to be the most reliable overall in estimating the wave amplitude and period. Specifically, it showed that, for the MSTID scenario, 80% of the AME is under 0.065 TECu (35% of the original amplitude) and for LSTID, AME is lower than 0.125 TECu (around one-third of the original amplitude). Among the techniques considered, FIF and SGOLAY also showed the highest reliability in preserving features in time; therefore, it is possible to conclude that FIF and SGOLAY are the best detrending/filtering techniques. In addition, if the computational time is a major constraint, SGOLAY would be the technique of choice since it is around 300 times faster than FIF.

  • Observational geometry (elevation and HIPP) can greatly impact the extracted wave period, especially for slowly moving TIDs. Specifically, we showed that the period error increases when lowering the elevation cut-off or increasing the ionospheric shell height (corresponding to the height at which the TID propagates in our 2D model). In the worst-case scenario, for MSTIDs with an elevation cut-off of 20° and HIPP equal to 350 km, period errors higher than 50% in 20% of the samples are expected.

  • The comparison between the reliable post-processing calibrated data and the NRT derivable observable has shown that the two verticalized observables, vTECNe and vTECPH assure a high degree of precision in the estimation of the amplitude. Under the LSTID scenario, the two vertical observables showed negligible AME (80% CDF lower than 0.01 TECu). In the MSTID case, instead, the NeQuick calibrated one (vTECNe) showed better performance, with an 80% CDF of 0.01 TECu against the 0.035 TECu for vTECPH.

  • Finally, this work affirms that a TID detection user can know beforehand the degree of amplitude and period error to expect. Moreover, we proposed and investigated some observables that can be leveraged for precise NRT TID detection and characterization, which, in our opinion, will soon be useful in correctly mitigating in real-time the TIDs disrupting effect on the growing number of GNSS and HF-based technologies.

Acknowledgments

The authors are grateful to Antonio Cicone (University of l’Aquila, Italy) for his valuable help with the FIF algorithm and useful discussions. In addition, we would like to thank Giorgio Savastano (Spire, Sapienza Università di Roma) and Mattia Crespi (Sapienza Università di Roma) for providing us with the VARION code and for the useful discussion related to the NRT observables and algorithms. The editor thanks Golekamang Thaganyana and an anonymous reviewer for their assistance in evaluating this paper.

Funding

Part of this work was supported by the Travelling Ionospheric Disturbances Forecasting System (T-FORS) project funded through the Horizon Europe programme of the European Commission under grant agreement 101081835, by the integrATed sysTEm for Multi-hazard from sPace over mediTerranean (ATTEMPT, https://progetti.ingv.it/it/pian-din#attempt-integrated-system-for-mutli-hazard-from-space-over-mediterranean) project (funded by the Ministry of University and Research – Fund for Investments of Central Administrations), and by the “Space Weather: PECASUS sistema di monitoraggio per applicazioni all “aviazione civile” project (funded by INGV).

Data availability statement

The Fast Iterative Filtering code is available at http://www.cicone.com. RINEX data from RING Network are available at http://ring.gm.ingv.it/. The arc database includes only 94 GPS satellites, and it is available at https://zenodo.org/doi/9510.5281/zenodo.10074542.

References

Cite this article as: Guerra M, Cesaroni C, Ravanelli M & Spogli L 2024. Travelling ionospheric disturbances detection: A statistical study of detrending techniques, induced period error and near real-time observables. J. Space Weather Space Clim. 14, 17. https://doi.org/10.1051/swsc/2024017.

All Tables

Table 1

Wave field parameters (amplitude, wavelength, azimuth, speed, and period) for the two different TID scenarios.

In the text

All Figures

thumbnail Figure 1

Reconstructed LSTID wave field calculated over Europe on a 0.5×0.5 lat/lon grid. The left panel represents the wave field calculated at t equal to 0, while the right panel represents the wave field at t equal to a half period. The black dot represents the location of the RING station AMUR.
In the text
thumbnail Figure 2

The top panels (a, b, c) show the different arcs of TEC considered in the reconstruction process under the MSTID scenario. Specifically, panel (a) shows the output of the GG calibration algorithm (vTEC), panel (b) shows the moving average smoothed vTEC arc and panel (c) shows the smoothed arc summed with the TID wave model (blue line of panel (d)). The bottom panel (d) instead shows the detrended (FIF-based band-pass) version of the synthetic arc (orange line) along with the TID wave model calculated at the IPP locations (blue line).
In the text
thumbnail Figure 3

The top panel (a) shows the different NRT-like observables studied together with the satellite elevation. The solid blue line is an example of a reconstructed arc, just like the one in Figure 2c. The orange line instead is the satellite elevation. The bottom panel (b) focuses instead on the detrended version of the four observables presented in panel (a), which share the same colour code. We zoomed into dTEC peaks corresponding to low, medium, and high elevations to better visualize the slight differences between the detrended observables.
In the text
thumbnail Figure 4

Amplitude error statistics for the six techniques considered. MSTID (right column) and LSTID (left column) parameters correspond to the ones presented in Table 1; the ionospheric thin shell height was set to 350 km. The top panels show the CDFs of AME for the different techniques (each line corresponds to the CDF of a given technique). X-axis shows the AME error in TECu, while the Y-axis corresponds to the probability P(AME) of having an error equal to or smaller than the corresponding X value. The bottom panels, instead, present the median along with their related percentiles (5, 16, 84, and 95) for the different filtering/detrending techniques.
In the text
thumbnail Figure 5

Violin plot for the correlation coefficient error. The X-axis shows the correlation error, while each violin (colour-coded according to different techniques) is a box plot where the shape represents the probability distribution of the correlation error for the given technique. The black diamond and line are, respectively, the median and mean of the correlation error distribution, while the shaded area represents the interquartile range. The left panel represents the LSTID scenario, while the right one is for MSTID.
In the text
thumbnail Figure 6

Cumulative distribution functions (CDF) of the Doppler error for different TID scenarios, heights of the IPP, and elevation masks. Plots on the right/left correspond to MSTID/LSTID, while top/bottom plots are for an IPP height of 250/350 km. The differently coloured lines show the CDF for different elevation masks, where blue, red, yellow, purple, and green are 20, 30, 40, 50, and 60 degrees, respectively. The X-axis represents the percentage period error, while the Y-axis shows the cumulative probability.
In the text
thumbnail Figure 7

Amplitude error statistics for the NRT observables derived dTEC. The top panels show the CDF of the amplitude error, with each line representing a different observable (blue, red, and yellow for sTECPH, vTECPH, vTECNe, respectively). The X-axis shows the amplitude error, while the Y-axis is the cumulative probability. The bottom panels show the median and the 18–84 and 5–95 percentiles for the different techniques. Here, the Y-axis represents the amplitude error in TECu. The plots on the left and right correspond to LSTID and MSTID, respectively.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.