‘ARTEMIS: Advanced Methodology Development for Real-Time Multi-Constellation (BDS, Galileo and GPS) Ionosphere Services’ Project Real-Time Ionospheric Services—Efficiency and Implementation

Abstract

This article discusses outcomes of the Polish–Chinese project dedicated to establishing multi-GNSS near-real-time ionospheric services. ARTEMIS (Advanced methodology development for Real-Time Multi-constellation (BDS, Galileo and GPS) Ionosphere Services) was a response to increasing GNSS data availability, including Galileo and BeiDou observations on the one hand and growing interest in high-quality ionospheric products on the other. The project resulted in elaborating methodologies to monitor the ionospheric Total Electron Content (TEC) and its fluctuations (ROTI index) based on a full multi-GNSS approach and establishing pilot real-time web services in a global and regional approach. The project’s outcomes are to be included in the International GNSS Service (IGS) and International Reference Ionosphere (IRI) in the near future. This article presents real-time ionospheric products developed under the ARTEMIS project and evaluates their performance using independent techniques such as DORIS observations and altimetry with regard to other existing products. The Discussion section also includes an evaluation of ARTEMIS products in positioning applications.

1. Introduction

The Earth’s ionosphere is an upper layer of the atmosphere, where cosmic mechanisms like EUV radiation and energetic particle precipitation stimulate ionisation of neutral atmospheric gases. The ionosphere starts around 50 km above ground and reaches up to the plasmasphere; however, its structure is highly diverse. As a result of ionisation and recombination mechanisms, ionospheric plasma has a characteristic vertical profile, including several regions with separate submaxima. Being magnetically susceptible, it also corresponds with the Earth’s magnetosphere and its characteristics; the ionosphere reveals serious dependency on the Earth’s magnetic field shape and strong latitudinal variability in the phenomena and mechanisms controlling plasma behavior. As ionisation and recombination balance is driven by solar energy, the ionosphere is also highly divergent in time—in very different scales. The temporal variability of plasma density includes diurnal, seasonal, and solar cycle-long periods. Furthermore, the variability resulting from magnetic field disturbances is also extremely important. The interplanetary magnetic field, magnetosphere, and ionosphere are connected in complex ways, and disturbances in the geomagnetic field often lead to disturbances in the ionosphere. The most rapid changes in the ionosphere occur during geomagnetic storms. Such events can significantly distort the daily changes in electron density in the ionosphere. Geomagnetic storms involve many different, interacting, and often opposing processes in the magnetosphere–ionosphere–thermosphere system.

Electromagnetic fields in the ionosphere are a serious threat for radio waves passing through. The Earth’s atmosphere blocks most of the electromagnetic spectrum, enabling lifeforms to exist on the planet’s surface. On the other hand, the selective transparency of the atmosphere limits the range of waves that can be utilised for the purposes of Earth–space communication, which nowadays plays an important role in telecommunication and navigation. Most of these links are based on radio waves. The ionosphere, due to its turbulent and dispersive properties for radio frequencies, is a source of the most serious and unpredictable errors in Earth–space communications.

A number of scientific entities, such as the International Reference Ionosphere (IRI) or International LOFAR Telescope (ILT), are interested in accurate and reliable information about the electron density within the ionosphere. One technology suffering from the impact of ionospheric plasma is global navigation satellite systems (GNSSs). Satellite positioning is based on measuring the distance between the receiver and at least four satellites simultaneously. Measurements are based on the L-band radio signals (around the level of 1.5 GHz). However, GNSS signals are disturbed by the ionosphere and this problem can be inverted. By using a network of permanent receivers with known positions, it is possible to retrieve information about the ionospheric delays of the signal and therefore about the ionospheric plasma density [1]. The International GNSS Service (IGS), for many years, has been providing global ionospheric maps (GIMs) [2]. With the growing availability of dense, world-wide permanent receiver networks with continuous operability, GNSSs became one of the main tools in ionosphere monitoring. For years, however, GNSS-derived products were mostly limited to the American Global Positioning System, NAVSTAR-GPS. Nowadays, there exist four operational GNSSs: the above-mentioned NAVSTAR-GPS (currently referred to as GPS), the Russian GLONASS, the European Galileo, and the Chinese BeiDou (often referred to as BDS). Most up-to-date receivers, including receivers installed within permanent networks (CORSs—Continuously Operating Reference Stations), support all four systems. It should be noted that simultaneous usage of all GNSSs requires addressing technical issues like intersystem bias (ISB) [3,4] and interchannel bias (ICB) in the case of the GLONASS system [5]. What is more, most permanent network receivers provide access to real-time data streams. The combination of four GNSSs significantly increases the number of available observations and, as a result, the quality of the resultant products. The development of real-time ionospheric products is strongly encouraged by both the Ionosphere and Real-Time committees of the IGS.

In response to a growing interest in ionosphere monitoring, not only in terms of a posteriori elaboration with several-day latency but also in near-real-time nowcasting approaches, the Space Radio-Diagnostics Research Centre (University of Warmia and Mazury in Olsztyn) together with the Academy of Opto-Electronics (Chinese Academy of Sciences) joined efforts to elaborate a near-real-time ionosphere monitoring service. The project ‘ARTEMIS: Advanced methodology development for Real-Time Multi-constellation (BDS, Galileo and GPS) Ionosphere Services’ was co-funded by the Polish National Research and Development Center and the Chinese Ministry of Science and Technology (MOST) within the first Polish–Chinese Joint Research Call.

Within the ARTEMIS project, Polish and Chinese teams developed real-time ionospheric products, including maps of ionospheric plasma total density. On the other hand, ARTEMIS services include maps of fluctuations of plasma density; many scientific parties are especially interested in ionospheric plasma turbulence and dynamics.

2. Materials and Methods

The ionosphere reveals dispersive properties in the radio frequency band; refractive index and therefore phase and group delays are dependent not only on the plasma density, but also on the signal frequency. The observation of at least dual-frequency signals then gives access to information about plasma density. All GNSSs are at least dual-frequency, supporting signals within bands like L1, L2, L5, etc. However, many simple navigation chips and receivers do not require multiple-frequency support; all professional receivers including those installed within permanent stations by default support all signals.

As mentioned above, ionospheric delay is present within GNSS signals. It is evident within so-called code (Equation (1)) and phase (Equation (2)) observation equations [6]:

$$P={\rho }_r^s+c(d{t}_r-d{t}^s)+T_r^s+I_r^s+c(d_r-d^s)+{ϵ}_r^s,$$

(1)

$$\Phi ={\rho }_r^s+c(d{t}_r-d{t}^s)+T_r^s-I_r^s+c({\delta }_r-{\delta }^s)+\lambda N+{ϵ}_r^s.$$

(2)

where P and $\Phi$ denote code and phase observations; $\rho$ is the geometrical distance between satellite s and receiver r; c stands for the speed of light; $dt$ represents clock offsets; T and I are tropospheric and ionospheric corrections; $d_r,d^s,{\delta }_r,{\delta }^s$ are instrumental delays in the receiver and satellite for code and phase observations, respectively; N is phase observation ambiguity; and $ϵ$ represents random errors. Some components such as tropospheric delays T and geometrical distance $\rho$ are frequency-independent and thus can be excluded from equations by a proper combination of dual-frequency observations—the so-called geometry-free (GF) linear combination:

$$P_{GF}=P_{L2}-P_{L1}.$$

(3)

$${\Phi }_{GF}={\Phi }_{L1}-{\Phi }_{L2},$$

(4)

where $P_{GF}$ and ${\Phi }_{GF}$ stand for geometry-free combinations of code (P) and phase observations ($\Phi$) with $L1$ and $L2$ carrier waves.

Ionospheric observables can be directly extracted from the geometry-free combination of GNSS dual-frequency pseudorange and carrier phase measurements. With the solution of Equations (3) and (4), besides the ionospheric component, there are also other frequency-dependent components: multipath, observation noise, and unknown integer phase ambiguities. To improve the quality of the derived ionospheric observables, the Carrier-to-Code Leveling (CCL) approach is used [7]. Within this approach, phase observations suffering from unknown but usually large integer ambiguities are used to smooth ambiguity-free but noisy code observations. Phase-smoothed code observations do not suffer from cycle slips, unknown ambiguities, and code noise. After solving other components (like including differential code biases (DCBs), remaining noise, etc.), the remaining ionospheric delay component holds plasma density as shown in Equation (5) [1]:

$$I=40.3{\displaystyle \frac{f_1^2-f_2^2}{f_1^2{f}_2^2}}{10}^{16}TEC,$$

(5)

where I is ionospheric delay (as in Equations (1) and (2)) and $f_1$ and $f_2$ are frequencies of $L1$ and $L2$. TEC stands for Total Electron Content. TEC is one of the main, general measures of ionospheric plasma density. It describes the column density of electrons—the amount of electrons integrated within a column with a 1 m² cross-section and a height equal to the signal’s path. TEC, mentioned in Equation (5), is then geometrically aligned with the satellite–receiver line of sight (LOS) and is often denoted as STEC (slant TEC). In order to address the TEC value independently of the satellite–receiver line of sight orientation, a relevant vertical value needs to be calculated. Currently, there exist many solutions with different approaches to the vertical plasma distribution problem; however, in most practical cases, the simple single thin-layer model of the ionosphere is adopted. This approach assumes the integration of the whole ionosphere into a single, infinitely thin shell at a certain altitude (usually 450 or 350 km). Within this approach, each LOS can be assigned to a single point at the ionospheric shell—the so-called Ionosphere Piercing Point (IPP). The dependency of observed slant and corresponding vertical TEC values is then described as follows [8]:

$$STEC=M\left(z\right)·VTEC$$

(6)

where

$$M\left(z\right)=[1-si{n}^{2}z·{(1+H_{ion}/R_E]}^{-1/2}$$

(7)

In Equations (6) and (7), z denotes the zenith angle of the LOS; $H_{ion}$—the altitude of the ionospheric thin shell; and $R_E$—Earth radius.

2.1. Ionospheric TEC Maps

As mentioned above, the key advantage of a GNSS in ionospheric monitoring is its global coverage with dense permanent receiver networks. For decades, GNSSs have been used as an efficient source of information for ionospheric map elaboration [2]. A number of authors still develop GNSS-based methods to monitor ionospheric activity at different scales in connection with events such as geomagnetic storms [9], but also in response to seismic activity [10,11]. The use of permanent receiver networks such as EUREF, CORS, or IGS gives an opportunity not only to monitor the ionosphere on a regional scale, but it also allows one to elaborate GIMs, which provide information about TEC in any place around the world. However, a serious issue is raised here. Permanent GNSS stations are inevitably constrained to physical land, whereas roughly 70% of the Earth’s surface is covered with seas and oceans. Lower-elevation observations can also provide information about electron content over the sea to a certain extent, but most of the ocean area is not covered by such measurements at all. In order to obtain a reliable, efficient global map, a proper interpolation algorithm needs to be incorporated. Such a technique should deal well with the problem of map ‘blank areas’ over the oceans, on the one hand, and fit properly to the local TEC variability described with receivers accumulated in dense networks, on the other hand. IGS-acclaimed product techniques include spherical harmonics expansion, Kalman filtering, and kriging [12].

Nowadays, many attempts being made to elaborate GIMs in real time. Efforts are guided by the IGS in the frame of its Ionospheric and Real-Time committees. IGS Ionosphere-Associated Analysis Centres (IAACs) are still developing real-time ionospheric products [13]; based on those, the IGS prepares a final, combined RT map. Currently, four centres provide RT products: Universitat Politècnica de Catalunya (UPC) [14], Centre National d’Etudes Spatiales (CNES) [15], Wuhan University (WHU) [16], and the Chinese Academy of Sciences (CAS). ARTEMIS products are developed under the auspices of CAS and UWM. CNES and WHU employ spherical harmonics to generate GIMs; UPC, until 2019, used kriging interpolation, but currently uses its own algorithm—Atomic Decomposition Interpolator of GIMs (ADIGIM) [13].

Real-time maps elaborated within the ARTEMIS project accommodate another technique used by IGS Ionosphere-Associated Analysis Centres (IAACs), namely the SHPTS method, which is used in GIM preparation by the Chinese Academy of Sciences (CAS). SHPTS stands for ‘Spherical Harmonic plus generalized Trigonometric Series function’. Ionospheric map generation with the SHPTS algorithm, described by Li et al. [17], can be divided into four steps:

• Based on the multi-GNSS observation from each permanent station, LOS STEC data are retrieved using phase smoothing of the code observations. It must be remembered, however, that after code observation phase smoothing, satellite and receiver differential code biases (DCBs) still remain within the equations and have to be properly handled.
• For each permanent station, a local ionospheric model is built. STEC values are calculated alongside corresponding vertical counterparts according to Equation (6). Local TEC variation is described using the Generalized Trigonometric Series function [18]. With this method, a two-dimensional (in terms of geographic latitude and local time) polynomial and a finite Fourier series (in terms of local time) are used to reconstruct the local TEC model. Within this model, receiver and satellite DCBs are also estimated with the IGGDCB method [19] (DCB solution method elaborated by the IGG (Institute of Geodesy and Geophysics), Chinese Academy of Sciences). The IGGDCB solution also addresses the problem of ISB.
• Global distribution of VTEC is then modelled with a spherical harmonic expansion (up to 15th order). Additional constraint on the spherical harmonic expansion is introduced to avoid obtaining negative TEC estimates in desolate locations over the oceans, where the inequality-constrained least squares (ICLS) method is applied. After deriving the initial weighted LS solution, negative regions are identified and non-negative constraints are imposed, ensuring physically realistic and reliable VTEC estimates [20].
• Combining the local and global models, the final GIM is created using a slightly modified DADS (Different Area Differential Stations) method [21]. This method fits local VTEC variability observed in dense observation coverage areas into the global model—in this case, built with spherical harmonics expansion, as described in step 3.

The SHPTS method is successfully used to produce RT-GIM (real-time GIM) products; however, the ARTEMIS project also involves regional elaborations. For the purpose of regional map elaboration, another method is incorporated, namely SHAKING (Spherical Harmonics Adding KrigING). Whilst SHPTS performs well in global elaborations, the SHAKING method is more sensitive for regional elaborations, where strong TEC gradients are expected—such as the equatorial zone.

In the SHAKING method, TEC modelling is split into two parts—deterministic and stochastic components. This is very important for elaboration in regions facing strong ionospheric gradients, as stochastic components become a dominant problem to solve. The regional deterministic part of TEC behavior is modelled with adjusted spherical harmonics (ASH), whereas the stochastic part is determined with the kriging method [22]:

$$I(\varphi ,\lambda )=ASH(\tilde{\varphi },\tilde{\lambda })+r(\varphi ,\lambda )$$

(8)

where

$$\left\{\begin{array}{c}{\displaystyle ASH(\tilde{\varphi },\tilde{\lambda })=\sum _{n=0}^{n_{max}}\sum _{m=0}^n{\tilde{P}}_{nm}sin\left(\tilde{\varphi }\right)·({\tilde{A}}_{nm}cos\left(m\tilde{\lambda }\right)+{\tilde{B}}_{nm}sin\left(m\tilde{\lambda }\right)}\hfill \\ {\displaystyle r(\varphi ,\lambda )=\sum _{i=1}^N{\lambda }_i{R}_{TEC\left(IP{P}_i\right)}}\hfill \end{array}\right.$$

(9)

where $I(\varphi ,\lambda )$ is the VTEC value depending on the geographic latitude $\varphi$ and longitude $\lambda$; $ASH(\tilde{\varphi },\tilde{\lambda })$ is the ASH-based deterministic component; and $r(\varphi ,\lambda )$ is the kriging-based stochastic component. ${\tilde{P}}_{nm}$ denotes the normalized associated Legendre function of degree n and order m; ${\tilde{A}}_{nm}$ and ${\tilde{B}}_{nm}$ are the ASH coefficients; N is the number of IPPs in a set range around ionospheric grid points; ${\lambda }_i$ is the corresponding weighting value; $R_{TEC\left(IP{P}_i\right)}$ is ionospheric residuals of IPP i. $\tilde{\varphi }$ and $\tilde{\lambda }$ are the adjusted latitude and cap longitude of the IPP [22].

2.2. TEC Fluctuation Index ROTI

The ionosphere reveals evident diurnal variability, as the dominant factor in ionospheric dynamics is the Sun. To accommodate this variability and guarantee a certain expected level of TEC continuity, linear interpolation in a 32 h time window is also included in the global modelling SHPTS algorithm. However, it should be remembered that the ionosphere is a highly turbulent medium. The number and complexity of factors cause ionospheric plasma to be unstable—especially in auroral, subauroral, and equatorial regions. Phenomena like energetic particle precipitation, plasma instabilities, acoustic gravity waves, and different coupling mechanisms trigger the generation of irregular structures within the plasma density. Such structures can cause unexpected behavior of radio signals passing the ionosphere. Ionospheric irregularities are often generated as a response to occasional, severe, hazardous events on the Sun or Earth. Technological and scientific communities are therefore interested not only in reliable, low-latency TEC information but also in plasma variability and dynamics description. Back in 2015, the IGS Ionosphere Working Group (currently Ionosphere Committee) introduced a new product describing TEC fluctuations—daily Rate of TEC Index (ROTI) maps for the Northern Hemisphere [23,24].

ROTI is a quantitative description of rapid TEC temporal fluctuations observed within a continuous satellite-tracking arc. Rate of TEC (ROT) is calculated as a time derivative of TEC series and expressed in TEC units (TECU) per unit of time (in most elaborations, minute intervals are adopted since such periods correspond with irregularities observed with GNSSs) [25,26]:

$$ROT={\displaystyle \frac{TEC\left(t_2\right)-TEC\left(t_1\right)}{t_2-t_1}}$$

(10)

ROTI is calculated as a standard deviation of ROT values collected within a selected period of time (in this manuscript, 5 min—similar to other studies) and spatial cell (in local elaborations, 1 degree by 1 degree in latitude and longitude) in order to describe the spatial distribution of TEC fluctuation occurrences:

$$ROTI=\sqrt{\langle RO{T}^2\rangle -{\langle ROT\rangle }^2}$$

(11)

3. Results

Efficient and reliable real-time ionospheric products are the main outcome of the ARTEMIS project. In the near future, these products will be fully available through two (Chinese and Polish) web services.

These include global TEC maps (GIMs) with an IGS-typical resolution of 2.5^∘ by 5^∘ in latitude and longitude and a time resolution of 5 min. Figure 1 presents an example of ARTEMIS GIMs; selected examples represent a whole day with a 2 h interval, but it should be remembered that ARTEMIS products are produced with 5 min intervals. The presented example describes ionospheric TEC behavior during an extreme geomagnetic storm that occurred on 10 May 2024.

In Figure 1, among other features, a serious deformation of the equatorial anomaly pattern is clearly visible in the evening. This is a well-known ionospheric behavior during high-geomagnetic-activity events. Such space weather events generate not only large-scale structures revealed by GIMs but also small-scale, rapid fluctuations in plasma density—especially at polar and subauroral latitudes. Those are successfully described with the ROTI index [27,28]. A number of authors [27,28,29,30,31] revealed ROTI enhancement and subauroral spread during geomagnetic storms. Therefore, real-time ROTI maps in global and regional scales were also prepared. Figure 2 presents regional European ROTI maps compared with corresponding TEC maps for the event described above. Figure 3 shows the evolution of disturbances in Earth’s magnetic field on 10 May 2024. Figure 4 presents corresponding European regional ROTI maps. Strong enhancement within ROTI values in both Figure 1 and Figure 4 can be clearly seen at high and even mid-latitudes (even below ${50}^{\circ }$N) after 18:00 UTC, when the geomagnetic storm began.

ARTEMIS also produced regional ionospheric maps using the SHAKING method, where the deterministic part is fitted by an ASH function and the stochastic component is generated by the kriging method. The SHAKING-based RT-RIM (real-time regional ionosphere map) demonstrates significantly higher accuracy compared to IDW and ASH methods, particularly during high geomagnetic activity, with improvements of up to 88% at edge stations [22]. This enhanced performance is attributed to SHAKING’s ability to correct distortions in edge areas and account for stochastic components, offering superior results across all geomagnetic conditions. Figure 5 presents an example of SHAKING RIMs above the Chinese sector.

ARTEMIS-elaborated real-time products are based on a multi-GNSS approach, since the number of observations is especially important with limited-observation processing. To better use the multi-GNSS observations and understand the effects of multi-GNSS data on ionospheric TEC modelling, we implemented different test scenarios, i.e., GPS + GLONASS, GPS + GLONASS + BDS, GPS + GLONASS + Galileo, and GPS + GLONASS + Galileo + BDS, for generating global ionospheric TEC maps. Figure 6 shows the accuracy of multi-GNSS generated global ionospheric maps compared to the IGS final GIM under different levels of geomagnetic activity. The results indicate that during periods of high geomagnetic activity, the fusion of multi-GNSS observation data does not significantly enhance the performance of GIM products. The model accuracy remains consistent across various configurations, likely due to a decrease in the precision of ionospheric observations from all systems during such periods. Conversely, during calm geomagnetic periods, models based on four-system observation data exhibit the highest accuracy, while solutions under the GPS and GLONASS dual-system mode perform the poorest. The accuracy of GIM during active geomagnetic periods is approximately 2–3 times lower than that during calm periods.

Similar to TEC GIMs, ROTI fluctuation maps were also elaborated using the multi-GNSS approach. In order to validate consistency between ROTI values obtained through observations from different systems, single-system (i.e., GPS only, GLONASS only, etc.) ROTI maps were produced. The 14 April 2022 strong geomagnetic storm was selected as a test period. Each elaboration (GLONASS, Galileo, BeiDou) was compared with the basic GPS product. Correlations between maps were obtained as follows: 0.7753 (GPS and GLONASS), 0.6776 (GPS and Galileo), and 0.7380 (GPS and BeiDou). Since no interpolation technique is involved in ROTI map elaboration, only ROTI cells co-occurring in both maps were taken for comparison.

4. Discussion

In order to assess and validate the real-time products’ reliability, a series of experiments were performed. Considering that reference data should have high accuracy and spatial coverage, empirical ionospheric models, such as the IRI model, were not selected as references. Within those, ionospheric information obtained from the different real-time products (those elaborated within the ARTEMIS project and others available) was compared with other ionospheric data sources:

• GIM VTEC—In the most straightforward approach, VTEC obtained from real-time products was compared with the IGS GIM—namely, the CODE rapid product CORG.
• DORIS dSTEC—Similarly to the GNSS dSTEC comparison, GIM-derived dSTEC was compared with observations made with Jason-3 DORIS. Thanks to the large relative frequency ratio between the two frequencies of DORIS, the theoretical precision of DORIS dSTEC is at the level of 0.028 TECu [32].
• GNSS dSTEC—Based on the real-time GIM VTEC, STEC LOS information was reconstructed and compared with dSTEC calculated from GNSS observations obtained from permanent receivers co-located with DORIS stations.
• Jason-3 VTEC—GIM VTEC was compared again with Jason-3 observation; however, in this approach, VTEC is extracted from altimeter data instead of DORIS. JASON-3 VTEC observations are considered reliable external references independent of GNSS and are frequently used in the assessment of GIMs, including IGS standards [33].

For each reference ionospheric data source, biases and standard deviations (STDs) have been calculated as follows:

$$\left\{\begin{array}{c}{\displaystyle Bia{s}_G=\sum _{i=1}^{N_{ref}}(TE{C}_{ref,i}-TE{C}_{G,i})/N_{ref}}\hfill \\ ST{D}_G=\sqrt{{\displaystyle \sum _{i=1}^{N_{ref}}{(TE{C}_{ref,i}-TE{C}_{G,i}-Bia{s}_G)}^2/(N_{ref}-1)}}\hfill \end{array}\right.$$

(12)

where $TE{C}_G$ is a proper TEC value (VTEC or STEC depending on the case) obtained from a GIM, $TE{C}_{r}ef$ is a reference value obtained from observations, and $N_{ref}$ stands for the number of reference samples.

The assessment comparison is first summarized in Figure 7 with the plot of GIM-TEC assessment bias and the STD results of CORG GIM versus time. A notable deviation can be observed, where RT-WHU (red line) and CNES (orange line) exhibit significantly worse performance than the RT-GIMs provided by the other three analysis centres. The STD of ARTEMIS, UPC, and IGS is at the level of 2.5–7.0 TECu, whereas it increases to 4.0–10.0 TECu for RT-GIMs of CNES and WHU.

The RT-GIM assessment process with DORIS STEC data was described in a previous work [32], where different RT-GIMs were cross-validated with data collected with 48 DORIS and 48 co-located GNSS receivers.

The process of calculating DORIS dSTEC closely resembles that of GNSS dSTEC. It relies on dual-frequency DORIS carrier phase measurements provided in RINEX DORIS format, analyzed along individual phase continuous arcs [32]. Both GNSS and DORIS dSTEC measurements carry ionospheric information, with GNSS reaching up to the Medium Earth Orbit (MEO) height and DORIS extending to Low Earth Orbit (LEO) heights. For the Jason-3 satellite, DORIS dSTEC observations cover altitudes of around 1300 km, encompassing the entire ionosphere and the majority of plasmaspheric electron contents (PECs). However, the influence of the remaining small portion of PECs is not factored into the assessment of DORIS dSTEC.

During DORIS STEC information retrieval, the DORIS antenna phase centre offset is also taken into account:

$$dSTE{C}_{DORIS}\left(t\right)={\displaystyle \frac{1}{40.3·(f_1^{-2}-f_2^{-2})*{10}^{16}}}\left[L_I\left(t\right)-L_I\left(t_{Emax}\right)-(\Delta D\left(t\right)-\Delta D\left(t_{Emax}\right))\right]$$

(13)

where $L_I$ stands for the geometry-free combination of $f_1$ and $f_2$ frequency observations, and $t_{Emax}$ denotes the epoch of the highest satellite elevation. $\Delta D$ is the associated geometry correction. The correction vector between the antenna mechanical reference point and the phase centre is [2.4128 m, −0.1325 m, 0.555 m] for the first frequency (2 GHz), and [2.4128 m, −0.1325 m, 0.9235 m] for the second frequency (400 MHz).

In order to evaluate the performance of different GIMs, each product (ARTEMIS, CNES, UPC, WHU, and final IGS) was compared to DORIS dSTEC observations. Daily values of mean bias and standard deviation (STD) values were used to compare products (Figure 8). The lowest daily STD values are recorded for the IGS product; values vary between 1.64 and 10.08. The UPC product reveals similar performance, with STD ranging from 1.51 to 10.49. For ARTEMIS, the WHU and CNES solutions’ daily STD varies between 2.10 and 11.32, between 2.30 and 12.17, and between 1.17 and 12.52. Average biases in regard to DORIS dSTEC for ARTEMIS, CNES, UPC, WHU, and IGS are 0.46, −0.21, 0.34, 0.56, and 0.22, respectively.

As mentioned above, a similar comparison was performed using 48 GNSS receivers co-located with DORIS stations [32]; again, biases and STDs are presented in the form of daily averages covering whole globe. The results are presented in Figure 9. Standard deviations of GIMs reveal similar performance of products. The average STDs of whole 3-year time series are 8.86, 11.04, 7.63, 12.70, and 8.66 for ARTEMIS, CNES, UPC, WHU, and IGS, respectively. The relative performance of particular GIMs is very similar to results obtained with DORIS. The UPC and IGS solutions also reveal the lowest biases: −4.06–4.07 with an average of −0.22 for UPC and −8.93–11.30 with −0.87 as the average for IGS. ARTEMIS RT-GIM biases range from −4.33 to 8.11 and the average bias is 0.90. ARTEMIS results are close to the IGS and UPC products’ performance. CNES reveals noticeable negative biases (with an average bias of −4.00 and values ranging from −14.64 to 10.51). This may suggest that the CNES product underestimates TEC values. Biases for WHU range between −7.67 and 13.30 with an average of 1.69. In 2023 and 2024, due to high solar activity bias, the STD levels are generally higher for all products.

The overall performance of individual GIMs versus DORIS and GNSS STEC observations confirms that the ARTEMIS product has similar accuracy to UPC and IGS combined maps, and slightly better than WHU and CNES.

Jason-3 VTEC has been proven to be a useful reference for VTEC measurements over the oceans. Jason-3 altimetry is equipped with a number of sensors that provide VTEC data independent of GNSS observations. As shown in the image, the orbit of Jason-3 can cover the entire ocean area from 65 degrees north to 65 degrees south latitude. The Poseidon-3B Altimeter operates at 13.575 GHz (Ku-band) and 5.3 GHz (C-band). Usually, we can obtain VTEC via the vertical phase ionospheric delay provided in the Ku-band frequency (Equation (14)). Since the orbit height of the Jason-3 satellite is about 1300 km, it can only measure the VTEC below 1300 km, resulting in a certain deviation between the VTEC information measured by Jason-3 and that measured by GNSSs. In this work, a sliding window of 20 s length was used to smooth the obtained VTEC time series to reduce the inherent noise effects.

$$VTE{C}_J=-{\displaystyle \frac{f_{Ku}^2}{40.3\times {10}^{16}}}{I}_{Ku}$$

(14)

where $I_{K}u$ is Ku-band ionospheric range correction and $f_{K}u$ is Ku-band frequency in GHz; $VTE{C}_{J,i}$ and $VTE{C}_{G,i}$ are VTEC extracted from Jason-3 and GIM observation i, respectively; $N_J$ is the number of involved observations.

The day-to-day bias and STD of each RT-GIM compared to Jason3 VTEC are calculated and plotted in Figure 10. Similar to the results generated in previous sections, UPC and IGS combined RT-GIMs show the best consistency, followed by ARTEMIS, CNES, and WHU. The STD of UPC and IGS products is about 3.0–6.0 TECu, which is 3.0–7.0 TECu for ARTEMIS, and 4.0–12.0 TECu for WHU/CNES, respectively. The RT-GIM from CNES shows negative deviations compared to other GIM products.

It can be seen that, although there are differences in detail under different evaluation methods, the average evaluation results are generally consistent. When applying these products, we recommend using a comprehensive approach of different evaluation methods to assess their accuracy. DORIS dSTEC and Jason VTEC, as two more independent methods, should be preferred.

4.1. RT ROTI Maps in Space Weather Event Monitoring

As mentioned before, ARTEMIS products are also part of space weather nowcasting systems. As Figure 11 shows, ROTI-described irregularity intensification and equatorward spread is well aligned with disturbances observed within the Earth’s magnetic field. SYM-H is a measure of longitudinally symmetric disturbances in the horizontal component of the magnetic field [34]. SYM-H is often treated as being equivalent to the Disturbed Storm Time (DST) index with a higher time resolution. Both SYM-H and DST are elaborated based on a chain of magnetometers at low latitudes; both indices describe the drop in the global magnetic field observed during geomagnetic storms.

ROTI observations at high and mid-latitudes reveal a great agreement with magnetometer-derived activity indices. ROTI nowcasting can therefore be an important tool in multi-instrumental space weather real-time monitoring. The limited number of available real-time data streams combined with multi-constellation tracking still gives satisfying resolution and efficiency in regional elaborations; ROTI intensification and equatorward spread is clearly visible in Figure 11.

4.2. RT VTEC Maps for Positioning Applications

The performance of the ARTEMIS real-time GIMs was also assessed in terms of satellite positioning impact. ARTEMIS product efficiency was compared with the BDGIM real-time product and IGS rapid (IGRG) maps. GIM performance was assessed in two positioning approaches: standard Point Positioning (SPP) and single-frequency precise point positioning (PPP). For the single-frequency standard point positioning (SPP), the ionospheric information obtained from GIMs is directly applied to correct ionospheric delay errors in GNSS pseudorange measurements (Equation (15)).

$$\left[\begin{array}{c}{P}_{r,f1}^1\\ P_{r,f1}^2\\ ⋮\\ P_{r,f1}^s\end{array}\right]-\left[\begin{array}{c}{I}_{GIM}^1\\ I_{GIM}^2\\ ⋮\\ I_{GIM}^s\end{array}\right]=B·\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\\ \Delta t_r\end{array}\right]$$

(15)

where $P_{r,f1}^s$ denotes the pseudorange observation equation of satellite s, frequency $f_1$ and receiver r; $I_{GIM}^s$ is the ionospheric correction obtained from the GIM for satellite s; B is the design matrix; and $\Delta x\Delta y\Delta z\Delta t$ correspond to position vector components. The equation actually represents the result after correcting satellite clock offsets, satellite biases, and tropospheric delay errors.

For single-frequency precise point positioning (PPP), the ionospheric product is used as the independent virtual observation and a certain prior variance is applied to constrain the real ionospheric estimates (Equation (16)).

$$\left[\begin{array}{c}{P}_{r,f1}^1\\ L_{r,f1}^1\\ P_{r,f1}^2\\ L_{r,f1}^2\\ ⋮\\ P_{r,f1}^s\\ L_{r,f1}^s\\ I_{GIM}^1\\ I_{GIM}^2\\ ⋮\\ I_{GIM}^s\end{array}\right]=H·\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\\ \Delta t_r\\ ZW{D}_r\\ I_{r,f1}^1\\ ⋮\\ I_{r,f1}^s\\ N_{r,f1}^1\\ ⋮\\ N_{r,f1}^s\end{array}\right]$$

(16)

where L corresponds to phase observations (similar to P in Equation (15)); H is the PPP design matrix; $ZW{D}_r$, $I_{r,f1}^s$, and $N_{r,f1}^s$ are network-based parameters for zenith wet delay (tropospheric component), ionospheric delay, and ambiguity parameter, respectively.

A total of 12 IGS reference stations with globally uniform distributions are selected to analyze the effects of different ionospheric model modifications on the performance of GPS L1 single-frequency SPP. The single-frequency SPP test covers four days, i.e., DOY 083, 084, 089, and 090. It should be pointed out that a strong geomagnetic event appeared on DOY 083 and 084, and the remaining two days experienced quiet geomagnetic conditions.

The RMS of residual positioning errors in the horizontal and vertical directions of single-frequency SPPS corrected by different ionospheric models is computed in comparison to the IGS precise coordinate solution. The ionospheric models for comparative analysis include ARTEMIS GIMs as well as the BDS-3 Global Broadcast Ionospheric Correction Model (BDGIM) [35] and the IGS rapid GIM product IGRG.

The overall performance of GPS L1 single-frequency SPP corrected by different ionospheric models is summarized in

Table 1. Single-frequency standard point positioning performance constricted with different ionospheric models (95% percentile).

	Ionospheric Model	Mean [m]	Min [m]	Max [m]
Horizontal	BDGIM	1.81	1.06	3.27
	ARTEMIS	1.22	0.69	2.06
	IGRG	1.19	0.46	2.65
	IGSG	1.08	0.40	2.01
Vertical	BDGIM	2.95	1.51	5.55
	ARTEMIS	1.87	0.75	3.81
	IGRG	1.90	0.79	3.36
	IGSG	1.65	0.60	3.10
3-D	BDGIM	3.46	1.84	6.44
	ARTEMIS	2.23	1.02	4.33
	IGRG	2.24	0.91	4.28
	IGSG	1.97	0.72	3.69

. For the horizontal component, the positioning accuracy of ARTEMIS RT-GIM and IGS rapid GIM products increased by 29.8% and 34.3%, respectively, compared to the BDGIM-corrected results. For the vertical component, the positioning accuracy of the above GIMs increases by 36.9% and 35.6%. The performance of RT-GIM is still slightly worse than the IGS rapid GIM generated in the post-processing mode. For a better understanding, results with the final IGS product—IGSG—were also included. Although the results are obviously better, ARTEMIS maps are not outperformed significantly. Mean errors are lower by 10% for IGS in regard to ARTEMIS.

Using the same stations for single-frequency SPP analysis, the quality of the generated RT-GIMs in the single-frequency PPP domain is analyzed in this section. Different from the direct ionospheric delay corrections in single-frequency SPP, the prior correction information provided by those ionospheric models is used to constrain the PPP model. In our analysis, the kinematic single-frequency PPP with additional ionospheric spatial and temporal constraints is applied. Similar to the previous single-frequency SPP analysis, the involved ionospheric models are ARTEMIS RT-GIM, BDGIM, and the IGS rapid GIM IGRG. The RMS of residual single-frequency PPP errors is used to validate the performance of different ionospheric models.

The IGS rapid GIM exhibits no significantly better performance than the RT-GIMs in the single-frequency PPP test. The overall performance of GPS L1 single-frequency PPP constrained by different ionospheric models is summarized in

Table 2. Single-frequency precise point positioning performance constricted with different ionospheric models (95% percentile).

	Ionospheric Model	Mean [m]	Min [m]	Max [m]
Horizontal	BDGIM	0.61	0.20	1.78
	ARTEMIS	0.46	0.12	1.42
	IGRG	0.48	0.15	1.35
Vertical	BDGIM	0.75	0.13	2.22
	ARTEMIS	0.62	0.16	1.75
	IGRG	0.60	0.16	1.64
3-D	BDGIM	0.97	0.24	2.85
	ARTEMIS	0.77	0.20	2.25
	IGRG	0.77	0.22	2.12

. For the horizontal component, the positioning accuracy of ARTEMIS RT-GIM and IGS rapid GIM products increased by 24.6% and 21.3%, respectively, compared to the BDGIM-corrected results. For the vertical component, the positioning accuracy of above three GIMs increases by 17.3% and 20.0%.

5. Conclusions

Nowadays, real-time ionospheric products play an important role both in positioning/navigation solutions and scientific purposes such as space weather monitoring.

Methods and products developed under the ‘ARTEMIS: Advanced methodology development for Real-Time Multi-constellation (BDS, Galileo and GPS) Ionosphere Services’ Polish–Chinese bilateral project present satisfying accuracy and usefulness for both scientific and application requirements.

The ARTEMIS products, while highly valuable in ionospheric monitoring and GNSS applications, have certain limitations that should be addressed. One potential area for improvement is the resolution and accuracy of the data, particularly in regions with sparse GNSS coverage, such as oceans or remote areas. To address this issue, the planned improvement involves incorporating DORIS and JASON observation data in future services. By integrating multi-source observational data, the GIM products will be updated to improve accuracy, especially over regions with sparse GNSS coverage, such as oceans.

ARTEMIS-developed products are also currently under IGS commissioning. All presented products will be available via dedicated WWW services. After appropriate assessment procedures, ARTEMIS products would also be introduced to IRI and ILT scientific communities.