1. Introduction
The ionosphere, which lies between 80 and 1000 km, plays an instrumental role in influencing the behavior of traversing electromagnetic waves due to its high electron density. The dynamic nature of the ionosphere leads to a variety of signal distortions, including attenuation, refraction, and delays, which are pivotal considerations for both terrestrial and extraterrestrial systems. Understanding these impacts requires a knowledge of the electron density distributions within the ionosphere, which is achievable through data inversion and reconstruction techniques. Various methods, including satellite observations, ionosondes, and incoherent radar systems have been employed to map the electron density distributions within the ionosphere [
1].
Computerized ionospheric tomography (CIT) is an effective tool developed to reconstruct ionospheric electron density values based on measurements that traverse the ionosphere [
2,
3]. The reconstruction techniques are typically based on iterative techniques or basis function expansions that represent the background ionosphere through the model matrix [
4]. Since it is not possible to take in situ measurements for the span of the 4-D ionosphere in latitude, longitude, height, and time, the measurements typically consist of satellite signals whose ray paths traverse the ionosphere, both in Low Earth Orbit (LEO) and in Middle Earth Orbit (MEO). The total electron content (TEC) is a pivotal parameter that is the line integral of electron density on a path as it corresponds to the total number of free electrons in a cylinder with a 1 m
2 cross section area [
5,
6]. The unit of TEC is TECU, where 1 TECU = 10
16 el/m
2. When TEC assessments are performed vertically, and aligning with the local zenith, the resultant metric is referred to as Vertical TEC (VTEC). Conversely, Slant TEC (STEC) denotes the TEC measurements that are conducted along trajectories deviating from the local zenith direction, and it typically encapsulates slanted paths. The STEC data become instrumental when processed through advanced tomographic calculations as they are designed to map the ionospheric electron density in three dimensions over time.
On the measurement side, the Global Navigation Satellite System (GNSS), notably Global Positioning System (GPS) (whose satellites are located at MEO), stands out as an important and cost-efficient instrument for ionospheric observations [
7]. By using the STEC estimated from the dual-frequency recordings between GPS satellites and both ground and LEO satellite receivers, valuable observations into the cumulative electron content along the line of sight can be achieved [
8]. Typically, radio occultation (RO) satellites at LEO, such as GPS/MET, CHAMP (Challenging Mini Satellite Payload), GRACE (Gravity Recovery and Climate Experiment), and FORMOSAT-3/COSMIC (Formosat Satellite 3/Constellation Observing System for Meteorology Ionosphere and Climate), have GPS receivers onboard [
8,
9,
10,
11]. The GPS-RO satellite paths enable one to trace the 3-D ionosphere in horizontal directions, whereas ground GPS receivers provide slanted paths in the vertical direction. The estimated STEC values from these combined paths alleviate the sparsity of measurement paths, enabling a detailed profiling of electron density across ionospheric layers.
Ionospheric tomography was first applied in 2-D using iterative methods such as Algebraic Reconstruction Technique (ART) [
2], Multiplicative Algebraic Reconstruction Technique (MART) [
12,
13], or Hybrid Algebraic Reconstruction Algorithm (HART) [
14]. While the advantages of these iterative techniques are their ease in implementation and low memory requirements, there are important disadvantages such as sensitivity to the level of noise, slow convergence, and dependency on the initial conditions.
In 3-D applications of tomographic reconstruction methods, one of the limiting factors is the scarcity of measurement data, resulting in ill-conditioned and sparse matrices in inversion [
15,
16,
17,
18]. Additionally, the inherent vertical bias of TEC measurements coupled with the uneven distribution of GPS receivers complicates the inversion process, leading to scenarios where some regions are oversampled while others lack sufficient data. Using ground-based observations alone, with their inadequate geometry for capturing vertical ionospheric structures, leads to models with limited resolution and accuracy. The addition of RO data complements ground-based measurements by enhancing vertical resolution and expanding coverage, especially in the upper layers of the ionosphere [
19,
20,
21,
22].
In addressing the challenge of data sparsity, various approaches have adopted a functional representation of electron density using a basis function expansion, and the problem has been transformed into estimating the coefficients for various bases [
4,
23,
24], such as Empirical Orthogonal Functions (EOFs) [
25], spherical cap harmonics (SCHs) [
26], Euclidean quadratic B-splines [
27], and Slepian functions [
28]. Achieving the desired spatial resolution across extensive areas using spherical harmonics, spherical cap harmonics, B-splines, and Slepian functions necessitates high numbers of coefficients to be estimated. This escalates memory usage and computational costs, thus presenting inversion challenges, particularly where GPS station coverage is sparse. Also, none of the basis representation studies have discussed the choice and calculation of the model matrix (which serve as a background for inversions), nor have they detailed the sampling matrix that was used for computation of the accurate path lengths.
IONOLAB-CIT is a 4-D regional electron density reconstruction method that integrates GPS-STEC data into a physical ionospheric model by optimizing the planar fit of the maximum critical plasma frequency of the F2 layer, foF2, maximum ionization height, hmF2, and surfaces in midlatitude ionosphere [
29,
30]. The method updates the ionospheric model by utilizing GPS-STEC, which is obtained from GPS recordings and fitting foF2 and hmF2 surfaces to a planar model utilizing only six parameters. The reconstructions are validated within a midlatitude region where the ionosphere exhibits a linear trend [
31]. To improve robustness, the outcomes of IONOLAB-CIT are temporally tracked and smoothed using Kalman filtering methods [
30]. IONOLAB-CIT provides the best electron density profiles in the midlatitude ionosphere; yet, for equatorial and polar ionospheres, where the trend structure deviates from planar, a more complicated model representation is required.
The singular value decomposition (SVD) of a model matrix provides one left and one right unitary matrix that contains the left and right basis vectors, as well as one diagonal singular value matrix whose singular values are in descending order, denoting the energy captured in the left and right singular bases [
32]. With a wide application area and a Digital Signal Processing approach, SVD can be used to capture the inherent physical structure of a model matrix by determining the total energy in each basis vector and singular value pair. SVD was applied to CIT for the first time in [
33] to obtain a basis decomposition of the global ionosphere. It has been shown that, when the model matrix is generated using a background ionospheric model and measurements using synthetic data, only four singular values are sufficient to represent the underlying physical structure of the ionosphere (even when there is additive noise in the measurement vector). For much smaller regions, SCH requires 48, B-splines requires 5832, and Slepian Basis Expansion requires 63 coefficients to be estimated [
26,
27,
28].
In this study, the proposed method in [
33] is further advanced to an SVD-based, closed-form CIT reconstruction algorithm for estimating the ionospheric electron density in 4-D, i.e., IONOLAB-Fusion. The generation of the model matrix is optimized by considering the solar activity and the hourly monthly variations in latitude, longitude, height, and time for any desired region. The sampling matrix is enlarged to include the STEC paths both between the ground receivers and GPS satellites, as well as between the receivers onboard the LEO satellites and GPS satellites. The developed algorithm also optimizes the STEC paths for uniform sampling in space by adding augmented virtual receivers in regions like oceans and deserts and decimating those in regions where there are too many in a close neighborhood to each other. The LEO satellite paths are computed within a temporal wide sense stationarity period and augmented for optimum reconstruction results. The measurement STEC can be obtained both from actual recorded GPS data or computed using the ionospheric model. The algorithm requires only one or two bases, and the reconstruction is completed in closed form. The algorithm was applied to Europe for the moderately solar active years for both quiet and disturbed days. When the reconstructed electron density was compared with those from the ionosondes and Global Ionospheric Map (GIM) TEC data in the region, it was observed that IONOLAB-Fusion has an excellent agreement with the ionosonde profiles and GIM-TEC when compared to those of the background ionospheric model.
In
Section 2, a brief methodological overview is provided for IONOLAB-Fusion.
Section 3 outlines the algorithmic application of IONOLAB-Fusion. In
Section 4, the results are presented, and then the conclusion is provided in
Section 5.
2. Computerized Ionospheric Tomography with IONOLAB-Fusion
The IONOLAB-Fusion algorithm, which is similar to the IONOLAB-RAY [
34,
35,
36,
37] and IONOLAB-CIT [
29,
30] algorithms, starts by dividing the ionosphere into 3-D spherical voxels and then sorting them lexicographically according to height, latitude, longitude and one epoch in time,
as illustrated in
Figure 1.
The
,
, and
in
Figure 1 are the starting points of the voxels, and
,
, and
are the resolutions for height, latitude, and longitude, respectively. Any [
] point in the grid can be defined as
where
,
, and
.
,
, and
denote the number of voxels in height, latitude, and longitude, respectively, which can be calculated as
where
is the initial hour,
is the final hour of the tomography, and
is the time resolution. The time index
can be defined as
where
, representing the index range of the number of time intervals included in one day. The
can be calculated as
Thus, in total, there are
entries considered in the reconstruction problem. The lexicographical index provides mathematical ease in handling the 4-D ionosphere in height, latitude, longitude, and time by numbering each voxel and ordering them into a 1-D vector [
3]. The lexicographical index
can be expressed as
The electron density defined for each voxel can be given as
where the subscripts
and
denote day and hour, respectively, and
is the transpose operator. To construct the model matrix and establish bases, the electron density calculated from an ionospheric model for
days is organized into a model matrix corresponding to the date and time of the tomography as
Singular value decomposition (SVD) is a powerful tool that is extensively utilized in CIT and mapping [
4,
14,
33,
38]. By decomposing the model matrix into physical basis functions, SVD can be used optimize the energy while minimizing the required number of physical basis vectors, thus making it suitable for the multivariate problems encountered in CIT. The SVD method decomposes
into three different matrices as
where
and
are the orthogonal matrices denoting the basis functions, and they are also the left and right singular matrices of Matrix
, respectively.
is a diagonal matrix, and it contains the singular values of the
as
where
represent the singular values in decreasing order as given below:
The singular values denote the energy in each corresponding basis vector in contribution to the model matrix. The energy content of the singular values in decreasing order allow for the separation of signal subspace and noise subspace. The singular values that correspond to the signal subspace can be identified by their significant contribution to the total energy in
. Thus, the term “significant” refers to the energy level of the singular values above a threshold. This threshold separates the signal subspace from the noise subspace. The number of significant singular values are denoted by
. By eliminating the matrix elements corresponding to the noise subspace, the signal can be de-noised. Thus,
to
in Equation (14) denote the signal subspace, and
to
denote the noise subspace. The electron density can be formed as the multiplication of the
with the appropriate coefficient vector
as
By examining the singular values,
, the ratio of the signal space energy to the total energy can be calculated as
By selecting the columns corresponding to the chosen singular values using
, the orthogonal basis functions of
can be determined. This ability to express high-energy phenomena with a minimal set of basis vectors enables SVD to effectively separate noise from useful information and to compress data, thus enhancing the accuracy and efficiency of ionospheric tomography. The significant electron density profiles corresponding to the signal subspace can be obtained by choosing only those significant singular values and basis vectors as
where
denotes the coefficients of the basis vectors.
The principle of CIT involves utilizing the measured STECs to represent the cumulative path lengths of the GPS signals passing through designated voxels, each of which are weighted by the electron density corresponding to the voxel.
is the number of satellite–receiver pairs at time
. Matrix
contains the distance in each voxel traveled by the rays between the GPS to ground-based receiver and the GPS to LEO-satellite-based pairs as they traverse through the voxels, and this can be referred as the sampling matrix. For the voxels that the rays do not pass through, zero is written in Matrix
as
where
is the length traversed in voxel number
, which lies in the path of the
th pair. The relationship between the STEC and electron density can be defined as
where
is the measurement vector that contain the STEC values. The measurement Matrix
can be defined as
The electron distribution can be expressed by basis components and reduces to the following form:
Coefficient estimation is performed in the least squares sense [
3,
4,
31,
33,
38] as
The estimated coefficients are then used together with the basis vectors to perform electron density reconstruction as
The application of the IONOLAB-Fusion algorithm is given in
Section 3.
4. Results
The IONOLAB-Fusion method, as outlined in
Section 3, was applied to Europe, a midlatitude region in the Northern Hemisphere, as an example to demonstrate the accuracy, robustness, and reliability of the algorithm. The European region consists of a dense ground-based GPS network and a significant number of ionosondes. Thus, the performance of the algorithm can easily be tested, and its validity can be established by using ionosonde electron density profiles and regional GPS-TEC maps. The region of interest spans latitudes between 40° N and 55° N and longitudes between −1° E and 15° E. The algorithm automatically expands the desired region by checking the possible slanted paths captured by the actual and virtual ground-based receivers. Thus, in the given scenario, the region of interest, which is identified with black dashed lines in
Figure 8, is expanded to encompass the STEC paths between the ground-based receiver stations situated within the region of interest and GPS satellites, as well as between the receivers onboard LEO satellites and GPS satellites. The boundary of the expanded region is situated between 34° N and 58° N latitude and −10° E and 25° E longitude, as shown in
Figure 8. This adjustment ensures that the oblique paths of the rays remain within the region throughout their trajectory, thus facilitating the accurate reconstruction of electron density. In
Figure 8, the ground-based receiver stations and ionosondes are denoted by red circles and blue triangles, respectively.
The ionosphere was partitioned into 3-D voxels with resolutions
and
in latitude and longitude, respectively. The range in height can be input by the user, employing variable or constant resolutions in km. The ionosphere is separated into distinct regions depending on the electron concentration, as discussed in detail in [
48,
49]. The maximum ionization takes place in the F2 layer, and the peak ionization height is denoted by hmF2. Typically, hmF2 takes values less than 400 km globally. The second important height boundary is the Chapman height, which is set to 428.8 km [
50,
51]. The ionization reduces above the Chapman height with a high gradient. While ionosondes reconstruct ionospheric density up to 1000 km, International Reference Ionosphere (IRI) provides electron density values up to 2000 km [
49,
52,
53]. The International Reference Ionosphere Extended to Plasmasphere (IRI-Plas) model can compute the electron density that exists up to the GPS satellite orbital height of 20,200 km [
40,
41,
42,
43,
44,
45]. Thus, depending on the choice of background ionospheric model, the user can define the range in height from 80 km to 20,200 km. In this study, we employed IRI-Plas-STEC to compute the contribution of TEC for various heights [
54]. IRI-Plas-STEC is provided as a space weather service at
www.ionolab.org (accessed on 10 November 2023). It was observed that, for the region of interest, the difference in the TEC between 2800 km and 20,200 km was less than 2 TECU. Therefore, we chose a height range of 90 km to 2800 km in this study.
The resolution in height can be selected as uniform in regular intervals or as non-uniform in correspondence with the significance of the electron density distribution. In this study,
km was used between 90 km and 600 km;
km was used between 600 km and 1300 km; and
km was used between 1300 km and 2800 km (similar to that in [
54]).
With the given spatial dimensions and resolutions, the ionosphere over the desired region is divided into spherical voxels, as discussed in
Section 2. The date and time of the reconstructions can be provided in a single entry to the algorithm. The temporal resolution has to be inputted in terms of hour. In this study, in order to keep the model matrix dimensions to a minimum,
was chosen to be 1. The lexicographical indexing helps to convert 4-D ionosphere into 1-D, as discussed in
Section 2 and shown in Equations (1)–(9). Finally, for
,
, and
we used
In this study, the formation of a model matrix was achieved through two distinct approaches. In the first approach, keeping the importance of solar radiation and solar activity on the ionosphere in mind, the model matrix was constructed based on an hourly monthly structure similar to the underlying ionospheric trend provided in the IRI and IRI-Plas models [
41,
44,
49,
52]. The solar activity can be grouped into four different categories depending on the Sun Spot Number (SSN) and solar flux (F10.7), which are the two main proxy indices describing the state of the sun [
55,
56,
57]. The years for which the annual average SSN was between 50 and 100, as well as annual average F10.7 that was between 100 and 130 sfu, were considered to be moderate activity years, as detailed in [
38] and the references therein. In the 23rd Solar Cycle (1996–2008), the 24th Solar Cycle (2008–2019), and the initial four years of the 25th Solar Cycle (2019–2030), 2003, 2004, 2011, 2012, 2013, 2015, and 2022 are classified as Moderate Solar Activity (MSA) years [
38]. In this study, 2003, 2004, and 2012 were chosen for the formation of the model matrix, whereas 2011, 2012, 2013, and 2015 were chosen for the test and validation years.
For the formation of an hourly monthly model matrix, we chose the month of April as an example. It was observed that the spring equinox season ended around the 10th to 14th of April, and the summer season started in the northern midlatitude ionosphere [
58,
59,
60]. Therefore, in the month of April, there are usually some severe geomagnetic storms, along with a large number of quiet days. In this study, for the first choice of model matrix formation,
was 90 days in total for the years 2003, 2004, and 2012.
The hours for reconstruction were chosen as 03:00 UT, 08:00 UT, 12:00 UT, 18:00 UT, and 23:00 UT, representing before sunrise, sunrise, noon, sunset and night hours, respectively. The model matrix was generated according to the procedure outlined in F1 and Algorithm 1 of
Section 3 for each hour separately as
.
In the second approach, the model matrix, , was formed by using the electron density obtained from the IRI-Plas model for the preceding days from the day of reconstruction. In this study, was set to 1 and 3 days.
The next step was the application of SVD to the model matrix
and
, as given in F2 and Algorithm 2 of
Section 3. The singular values
and the normalized cumulative energy
for
were formed by three MSA years and the days in the month of April. These are provided in
Figure 9, up to
10 for all hours of investigation, as an example. As shown in
Figure 10a, it is evident that the major contribution was due to the first singular value representing the signal subspace for all hours. The first singular values for all hours exceeded
, whereas the second singular values had values smaller than
. The energy of the first basis ranged between 97.5% and 99% across all hours, covering a significant amount of the signal subspace. Since the difference between the energy of the first and second bases was minimal, indicating a small second singular value compared to the first one, the contribution of the second basis was negligible. Thus, for all hours that were considered, only one significant basis was enough to capture the trend. For the disturbed days during sunset and sunrise hours where the ionospheric variability increases, the basis number can be increased to two depending on the location of the ionosonde in comparison.
was set to 1 and 3 days for the formation of , and the significant singular values were determined to be one for and two for
In this study, the example days for the tomography were chosen as 6 April 2011, 17 April 2011, 24 April 2012, 4 April 2013, 25 April 2013, and 8 April 2015. The geomagnetic indices of the Planetary K-Index (Kp), Disturbance Storm Time (Dst) (nT), and Auroral Electrojet (AE) (nT) for the abovementioned dates are provided in
Figure 10a, 10b, and 10c, respectively. The data for the plots were obtained from
https://omniweb.gsfc.nasa.gov/form/dx1.html (accessed on 18 October 2023). The day boundaries are given in vertical red dashed lines. The disturbance indicators are provided in horizontal black dash dotted lines. The Planetary K-Index (Kp) index serves as a worldwide gauge, indicating the level of geomagnetic activity within a specific three-hour timeframe [
44,
45,
61]. The Kp values larger than 3 are indicators of a geomagnetic disturbance. The Disturbance Storm Time (Dst) index is a measure of the magnetic storm level [
44,
45,
61]. It is computed from the horizontal variations of the geomagnetic field that is measured at four low-latitude stations that are distributed globally in longitude. When the Dst is smaller than −30 nT, a geomagnetic disturbance is considered as set. The Auroral Electrojet (AE) index is derived from the geomagnetic variations in the horizontal component of the geomagnetic field that exists along the auroral zone in the Northern Hemisphere [
44,
45,
61]. A day is considered to be disturbed when the AE is greater than 500 nT [
61]. According to these disturbance indicators, 17 April 2011, 4 April 2013, and 8 April 2015 are geomagnetically quiet days, and 6 April 2011, 24 April 2012, and 25 April 2013 are disturbed days with geomagnetic storms.
In the densely populated ground-based receiver region of Europe, we conducted a study to prevent overlapping rays between the receivers and GPS satellites within the same voxels, thus mitigating the risk of overfitting. This was achieved by selectively reducing the number of stations, ensuring a uniform distribution after decimation. In the IONOLAB-Fusion algorithm, ground-based receiver stations in the region of interest were detected automatically, and a homogenously distributed subset was identified with a rate of about 20% using uniform sampling. Consequently, for the example region of this study, 48 almost uniformly distributed ground-based receivers out of possible 163 GPS receivers were chosen, as presented in
Figure 8 with red circles.
The GPS satellites visible to these receivers at elevation angles exceeding 40° are identified automatically.
Figure 11 illustrates the rays between the ground-based receiver station and GPS satellite pairs for 8 April 2015 at 12:00 UT. In
Figure 11a, 3-D representations of rays between pairs are provided. The selected region is indicated as the patch and the GPS satellites are identified by their Pseudorandom Noise (PRN) numbers.
Figure 11b displays the projection of the 3-D representation onto a 2-D latitude and longitude plane.
Figure 11c depicts a zoomed-in plot focused on the boundary of the extended region.
The COSMIC-1 mission, launched on 15 April 2006, stands as the world’s first multi-satellite occultation mission [
62]. Comprising six identical satellites, each positioned at 72° and an altitude of 810 km, COSMIC-1 is equipped with advanced GPS radio occultation receivers onboard. In this study, the position of the LEO satellites from the COSMIC-1 mission was used. If the user opts to utilize the rays between the receivers onboard the LEO satellites, the GPS satellites within the line of sight of the LEO satellites are identified automatically.
In this study, the optimal temporal resolution between the LEO and GPS satellite pairs was obtained by considering the wide sense stationarity (WSS) period of the midlatitude ionosphere, whose details are given in [
63,
64,
65,
66]. WSS is the duration where the first and second moments of the experimental probability density function remain unchanged. In the studies of [
64,
65,
66], the WSS period of the midlatitude ionosphere was determined to be 11 min. Within the WSS period, the temporal resolution is chosen by taking two conflicting constraints into consideration. It is desirable to take as many samples as possible within the WSS period, and, at the same time, the slanted ray paths should not cross the same voxels for the given spatial resolutions in the formation of 3-D spherical voxels. Thus, a temporal resolution of 30 s was selected to prevent voxel overlap, thus ensuring that the rays from different pairs traverse distinct voxels. At the same time, 22 samples can be obtained within 11 min WSS periods. In
Figure 12a, an example of the 3-D representations of the rays between GPS and LEO satellites for 8 April 2015, at 12:00 UT, during a 11 min WSS period with 30 s temporal resolution is given.
Figure 12b displays the projection of the 3-D representation onto a 2-D latitude and longitude plane.
Figure 12c depicts a zoomed-in plot focused on the boundary of the extended region.
The STEC data between the receiver and GPS satellite pairs were calculated to form the measurement vector in Equation (19). In this study, the GPS-STEC values were estimated using the IONOLAB-STEC algorithm given in [
5,
6,
51] with a 30 s time resolution. The receiver Differential Code Biases (DCBs) were estimated using the IONOLAB-BIAS algorithm provided in [
67]. IONOLAB-STEC was conveniently accessible as an online space weather service via
www.ionolab.org (accessed on 10 November 2023) [
68].
The STEC between the LEO and GPS satellites was computed using the IRI-Plas-STEC algorithm [
54]. Online computation of the IRI-Plas model is provided at
www.ionolab.org (accessed on 13 November 2023) [
69]. In the studies of [
70,
71], it is shown that, for ground-based receiver station and GPS satellite pairs, IRI-Plas-STEC values are very similar to IONOLAB-STEC values for quiet days in the midlatitude ionosphere.
The IONOLAB-Fusion algorithm allows the user to define virtual slanted ray paths that traverse the ionosphere in a horizontal direction, as outlined in Flowchart F3.a and Algorithm 3.a of
Section 3. For the example scenario given in
Figure 11 and
Figure 12 on 8 April 2015 at 12:00 UT, the number of actual ground to GPS satellite IONOLAB-STEC ray paths was 172; the number of actual GPS-LEO IRI-Plas-STEC ray paths in a WSS period was 41; and the number of virtual augmentation ray paths was 25.
In this study, the paths of the rays between ground-based and LEO-based receivers and GPS satellites were traced as they traversed the ionosphere using the wave propagation algorithm IONOLAB-RAY to generate the sampling matrix in Equation (18). The IONOLAB-RAY algorithm uses the IRI-Plas model to obtain the physical parameters of the ionosphere, such as the electron density. The sampling matrix is formed by calculating the path lengths in the voxels for the pairs, as given in F3.b and Algorithm 4 of
Section 3. If desired, the IONOLAB-RAY algorithm can automatically input Global Ionospheric Map (GIM)-TEC values to update the background IRI-Plas model to the current state of the ionosphere. The GIM-TEC maps with a temporal resolution of 2 h and a spatial resolution with 2.5° in latitude and 5° in longitude are automatically interpolated in the IONOLAB-RAY algorithm for the user-selected resolutions in terms of latitude, longitude, and time. In this study, the GIM-TEC maps were obtained from the Center for Orbit Determination in Europe (CODE) (
ftp://cddis.gsfc.nasa.gov/pub/gps/products/ionex/ (accessed on 13 November 2023)).
As depicted in F4 and Algorithm 5 of
Section 3, the measurement matrix was formed using the sampling matrix and subspace basis vectors. The association between the measurement vector and the measurement matrix was established using basis coefficients. The basis coefficient vector was estimated in the least squares sense, as shown in Equation (22). Ionospheric electron density reconstruction was estimated using the basis coefficient vector and signal subspace basis vectors for
voxels for each time epoch in a closed form. In
Figure 13, we present slice representations of the electron density reconstructions within a selected region of the ionosphere for both quiet and disturbed days at different times, as displayed in the leftmost column. The middle column shows the electron densities that were calculated using the IRI-Plas model for the same dates and times. The rightmost column depicts the differences between the electron density reconstructions and the densities calculated from the IRI-Plas model. The reconstruction results utilized only the ground-based receivers in the IONOLAB-Fusion algorithm. As shown in
Figure 13a,d, the electron density was reconstructed using Model Matrix
for the ionospheric disturbed day of 25 April 2013, 03:00 UT and the ionospheric quiet day of 17 April 2011, 08:00 UT, respectively. As shown in
Figure 13g, Model Matrix
was formed using
1 for the ionospheric quiet day of 8 April 2015, 18:00 UT. As shown in
Figure 13j, Model Matrix
was formed using
3 for the ionospheric disturbed day of 24 April 2012, 23:00 UT. The slice heights were selected between 100 km and 500 km with a resolution of 100 km.
In
Figure 13a, the electron density is shown for 25 April 2013 at 03:00 UT. During night hours, the ionization levels were low. The highest electron density was observed at 300 km in the southeast part of the region, followed by 400 km in the southwest part of the region.
Figure 13d displays the electron density for 17 April 2011 at 08:00 UT. As the sun rises, the ionosphere begins to ionize, leading to an increase in electron density. The eastern part of the region exhibits higher electron density values due to the sunrise. The highest electron density is seen at 300 km on the east side of the region, with the lower altitudes also starting to ionize. In
Figure 13g, the electron density for 8 April 2015 at 18:00 UT is provided. Compared to 08:00 UT, the electron density at 18:00 UT remained higher due to the cumulative effect of the solar radiation throughout the day. The highest electron density was at 300 km, spreading across the region. In
Figure 13j, the electron density is shown for 24 April 2012 at 23:00 UT. Since ionization depends on solar radiation, after sunset, the electron density decreases at all altitudes. As the sun sets in the West, electron density values are higher in the western part of the region. The higher electron density depicted in
Figure 13a when compared to
Figure 13j can be attributed to the presence of a storm during the early hours of the day. Since there was a disturbance in the early hours of 25 April 2013, the electron density of that day, as reconstructed in
Figure 13a, captured the increase in the disturbance more accurately than the electron density calculated using the IRI-Plas model shown in
Figure 13b, as indicated by the positive difference in
Figure 13c. On quiet days, as depicted between
Figure 13d,i, the electron densities captured by the IONOLAB-Fusion and IRI-Plas approaches showed the most significant differences at 300 and 400 km slices, which is where the peak electron density occurs, as given in
Figure 13f,i. There was a depletion in the electron density on 24 April 2012 at 23:00 UT due to a negative disturbance, which was detected by IONOLAB-Fusion, and this resulted in a negative difference that is visible in
Figure 13l.
To validate the results obtained using IONOLAB-Fusion, the electron density values at the location of the ionosondes were compared with the ionosonde vertical electron density profiles. The ionosondes used in this study, along with their respective locations, are detailed in
Table 1.
The normalized metric distance was obtained by calculating the ratio of the metric distances between two vectors to the reference vector as
where
and
are two vectors. Symmetric Kullback–Leibler Distance (SKLD) provides the entropical distance between two probability density functions and compares the formal similarity [
72]. The SKLD is calculated as
where
, the empirical probability density function, is defined as
The height at which the ionization reaches its peak is denoted as hmF2, while the corresponding peak value is defined as NmF2 [
48]. The height at which the Single-Layer Ionosphere Model (SLIM) assumes that the ionosphere is an infinitely thin layer is the Chapman height, which is considered to be at 428.8 km [
50,
51]. Since the ionosonde provides the reconstructed vertical electron density profile values up to 1000 km using the near vertical incidence sounding ionograms, comparisons were determined up to 1000 km. The examples of the electron density profiles at the ionosonde positions were calculated from the electron density reconstructions given in
Figure 13, and they were then plotted alongside the vertical electron density profiles obtained from the ionosondes and the vertical electron density profiles calculated by the IRI-Plas model, as shown in
Figure 14. In
Figure 14, the following electron density profiles are given: the ionospheric disturbed day of 25 April 2013 at 03:00 UT at the position of DB049 in
Figure 14a; the ionospheric quiet day of 17 April 2011 at 08:00 UT at the position of PQ052 in
Figure 14b; the ionospheric quiet day of 8 April 2015 at 18:00 UT at the position of JR055 in
Figure 14c; and the ionospheric disturbed day of 24 April 2012 at 23:00 UT at the position of EB040 in
Figure 14d.
As shown in
Figure 14, the electron density that was reconstructed with IONOLAB-Fusion using the IRI-Plas model aligned better with the ionosonde vertical profiles up to hmF2. However, at higher altitudes, which is known as the topside, differences in the underlying models between the IONOLAB-Fusion and ionosonde inversion algorithms led to discrepancies between the reconstructed electron density values and ionosonde vertical electron density profiles [
73]. This was especially the case beyond the Chapman height, where notable disparities emerged between the ionosonde vertical electron profiles and the reconstructed electron density profiles that were constructed by IRI-Plas.
The performance of the electron density reconstruction by IONOLAB-Fusion was evaluated by comparing the reconstructed electron density values at the ionosonde positions with the vertical electron density profiles of the ionosondes, and this was achieved by employing the comparison methods on different heights. The comparison methods were separately applied to the hmF2, Chapman, and 1000 km heights. Additionally, to assess the improvement of the tomography results compared with the ionospheric model, the electron density calculated by the IRI-Plas model at the ionosonde positions was also compared with the ionosonde vertical electron density profiles as reference values. The comparison results of
Figure 14 are given in
Table 2.
On the ionospheric quiet day of 17 April 2011 at 08:00 UT, at PQ052, the Model Matrix
was utilized and the electron density was reconstructed using IONOLAB-Fusion. It was observed that the reconstructed electron density closely resembled the vertical electron density of the PQ052 across all altitudes when compared to the IRI-Plas model. The similarity in SKLD values between the electron density reconstructed by IONOLAB-Fusion and the electron density calculated by the IRI-Plas model suggests that both models exhibited a similar shape to the vertical electron density of the PQ052 up to the Chapman height. However, both profiles differed from the vertical electron density profiles of the ionosonde after the Chapman height. Despite this, at 1000 km, the reconstructed electron density outperformed the one calculated by the IRI-Plas model. On the quiet day of 8 April 2015, at 18:00 UT, the model matrix
is formed by using
1. Similarly, the reconstructed electron density at the hmF2 and Chapman heights closely matched the vertical electron density profile of the JR055 when compared to the electron density calculated by the IRI-Plas model. However, due to the model matrix only containing a priori information from the previous day, the reconstructed electron density was unable to accurately track the ionosonde vertical profile after hmF2. On the disturbed day of 25 April 2013 at 03:00 UT, Model Matrix
was utilized at DB049. The reconstructed electron density profiles exhibited a closer resemblance to the DB049 vertical electron profile in both value and shape compared to the electron density calculated by the IRI-Plas model at all altitudes, as indicated by the results with smaller values obtained from both the metric distance and SKLD methods. On the day of 24 April 2012 at 23:00 UT, at EB040, the ionosonde vertical electron density values utilized for normalization were relatively smaller compared to those depicted in
Figure 14a,b. Consequently, the minor changes had a more pronounced impact, leading to larger metric distances. Moreover, the inability of the ionospheric models to precisely depict the ionosphere during disturbed days resulted in larger disparities in comparison to the method results for electron density profiles derived from both the IONOLAB-Fusion method and the IRI-Plas model. On the disturbed day of 24 April 2012 at 23:00 UT, Model Matrix
was formed using
3. The reconstructed electron density profile better reflected the peak electron density value, NmF2, compared to the IRI-Plas model. However, due to the disturbance, the prior information in the model matrix could not accurately represent the current ionosphere. Consequently, the normalized metric distance and SKLD values were higher at this time. While the reconstructed electron density yielded better results in terms of metric distance, the SKLD values were larger due to the shift in the hmF2 height.
To calculate the overall results, Model Matrix
was formed using
days for the MSA years for April, while
was formed using
3 previous days for the specified dates and times. Electron density reconstruction was conducted by IONOLAB-Fusion utilizing both
and
for all the specified dates and times. The reconstruction involves ground-based receiver station and GPS satellite pairs, as well as LEO-GPS satellite pairs and virtual LEO paths.
Table 3 presents a comprehensive comparison of the tomography results and all the ionosonde vertical electron density profiles, which encompass all the quiet and disturbed days across all hours. Furthermore, the electron density values computed from the IRI-Plas model were compared with the same dataset, and they were also used as a reference. To assess the impact of the horizontal paths, reconstructions were performed with and without them in various scenarios. The scenario column in
Table 3 indicates the dataset utilized in the IONOLAB-Fusion method for electron density reconstruction. In Scenario 1 (S1), reconstruction was performed using paths between only the ground-based receiver stations and GPS satellite pairs. In Scenario 2 (S2), the horizontal ray paths between the LEO-GPS satellite pairs at the epoch of concern were added to the data from S1. In Scenario 3 (S3), the paths between the LEO-GPS satellite pairs within the WSS period were added to the data from S1. Finally, in Scenario 4 (S4), virtual LEO rays were added to the data from S3, which comprise both horizontal and slanted actual, virtual, and augmented ray paths.
In S1, upon analyzing the comparison results up to hmF2 for quiet days, it was evident that the reconstruction results that used and exhibited greater similarity with the vertical electron density profile of the ionosonde compared to the electron density values calculated from the IRI-Plas model. Moreover, reconstructions utilizing as the model matrix outperformed those using in terms of the normalized metric distance. Examination of the SKLD results revealed that the IRI-Plas model results exhibited a similar shape to the ionosonde vertical electron density values compared to the reconstruction results utilizing . This similarity arises from the IRI-Plas model’s consistent performance in the midlatitude region during quiet days. Conversely, the reconstruction results utilizing yielded SKLD values closer to those of the IRI-Plas model. This can be attributed to the incorporation of more a priori information in the model matrix. The comparison results of and were identical for the quiet days across all scenarios. In S2, there is a slight enhancement in the results of the , , , and up to the hmF2 height with the addition of momentary RO data; however, this improvement was not substantial. Furthermore, incorporating the RO data within the 11 min WSS period led to an enhancement in the comparison method results of and . While there were improvements in the and results in the S3 scenario, which utilized virtual LEO rays, there were also improvements in , , , and up to the hmF2 height. When analyzing the comparison of results up to the Chapman height and 1000 km for quiet days, it was noted that the electron densities that were reconstructed utilizing and produced better results in metric distances up to the Chapman height across all scenarios. However, at 1000 km, larger comparison results were obtained, specifically ones that are similar the results obtained by IRI-Plas. This discrepancy arises from the increased disparities between the ionospheric model utilized by IONOLAB-Fusion for reconstruction and the ionosondes employed for the inversion algorithm. These results suggest that, on quiet days, employing Model Matrix generates electron density profiles, even without the use of RO data. Moreover, without the use of RO data, Model Matrix , which integrates more prior information, outperformed the results obtained from . Additionally, the incorporation of the RO data obtained from within a WSS period and virtual LEO rays enhanced the metric distance results up to the Chapman height for the tomography using both and , with more pronounced improvements being observed when utilizing . Upon analysis of the results for the disturbed days up to the hmF2 height, Chapman height, and 1000 km, it was evident that electron density reconstruction using both and outperforms the electron density calculated by the IRI-Plas model. In particular, the results obtained from Model Matrix , which integrates more prior information, exhibited a superior performance when compared to . Scenarios S2 and S3 were not considered as there was no LEO ray pass through the region on disturbed days.
The TEC values for all the latitude and longitude points within the region could be calculated by multiplying the electron density within each voxel by the corresponding increment of the voxel height and then summing these products over the entire height range. The routinely provided GIMs from the IGS are widely regarded as a dependable ionospheric product. In this study, the GIMs sourced from the Center for Orbit Determination in Europe (CODE) were employed to assess the quality of the ionospheric TEC calculated by IONOLAB-Fusion. Additionally, the TEC calculated from the IRI-Plas model was compared with GIM-TEC to serve as the reference value. The TEC was computed from the IONOLAB-Fusion electron density profiles using the IRI-Plas-STEC algorithm [
68] for each grid point in latitude and longitude for the reconstruction height. In
Figure 15, examples of the TEC maps from CODE-GIM, IONOLAB-Fusion-TEC, and IRI-Plas-TEC are presented for the ionospheric quiet day of 8 April 2015 at 08:00 UT. CODE-GIM is typically provided with resolutions of 2.5° in latitude and 5° in longitude, and it is interpolated to a 1° resolution in both the latitude and longitude of the reconstruction. The maps are provided for the extended region that encompassed the region of interest, as shown in
Figure 8.
As shown in
Figure 15, the TEC map derived by IONOLAB-Fusion exhibited closer similarity to the GIM-TEC results both in terms of value and the distribution of electron density when compared to the TEC map generated from the IRI-Plas model. The difference between the maps can be calculated using the root mean square error (RMSE) as
where
denotes the GIM-TEC and
represents the calculated TEC maps.
The RMS values were computed for all quiet and disturbed days of reconstruction, through all hours, and for all the grid points between CODE-TEC and the maps calculated by the IONOLAB-Fusion and IRI-Plas models. The comparison results are given in
Table 4. For the extended region in
Figure 8, the
and
values in Equation (27) were 25 and 36, respectively, making 900 grid points in total. For the height, electron density computed up to 2800 km was used in the IRI-Plas-TEC and IONOLAB-Fusion-TEC maps.
It was observed that, when the GIM-TEC was taken as the ground truth, the TEC maps calculated with IONOLAB-Fusion performed excellently compared to those obtained from the IRI-Plas model, with improvements of 37.89% for quiet days and 31.58% for disturbed days (even though the full ionospheric and plasmaspheric heights up to the GPS satellite were not used).
The IONOLAB-Fusion algorithm, with its ease in implementation, accuracy in comparison with ionosonde and GIM-TEC, and reliability and robustness in application to any region on Earth, provides a versatile means of estimation and a method through which to update electron density profiles for any ionospheric condition. It was observed that augmentation of virtual rays within the wide sense stationarity period of the ionosphere enhanced the accuracy of the reconstruction results when compared to ionosonde and GIM-TEC. The overall excellent performance is evidence that IONOLAB-Fusion is a strong candidate for an effective method to reconstruct electron density profiles for any date, time, region, and ionospheric state.