An Extended Simultaneous Algebraic Reconstruction Technique for Imaging the Ionosphere Using GNSS Data and Its Preliminary Results

: To generate high-quality reconstructions of ionospheric electron density (IED), we propose an extended simultaneous algebraic reconstruction technique (ESART). The ESART method distributes the discrepancy between the actual GNSS TEC and the calculated TEC among the ray–voxels based on the contribution of voxels to GNSS TEC, rather than the ratio of the length of ray–voxel intersection to the sum of the lengths of all ray–voxel intersections, as is adopted by conventional methods. The feasibility of the ESART method for reconstructing the IED under different levels of geomagnetic activities is addressed. Additionally, a preliminary experiment is performed using the reconstructed IED proﬁles and comparing them with ionosonde measurements, which provide direct observations of electron density. The root mean square errors (RMSE) and absolute errors of the ESART method, the simultaneous algebraic reconstruction technique (SART) method, and the International Reference Ionosphere (IRI) 2016 model are calculated to evaluate the effectiveness of the proposed method. Compared to the conventional SART method of ionospheric tomography and the IRI-2016 model, the reconstructed IED proﬁles obtained using the ESART method are in better agreement with the electron density obtained from the ionosondes, especially for the peak electron densities (NmF2). In addition, a case study of an intense geomagnetic storm on 17–19 March 2015 shows that the spatial and temporal features of storm-related ionospheric disturbances can be more clearly depicted using the ESART method than with the SART method.

The unstable ill-posed problem is crucial for GNSS-based CIT techniques due to the sparsity of GNSS tracking stations and the limited viewing angles from GNSS observations [23][24][25][26], which leads to difficulties in accurately constructing the IED distributions for GNSS tomography. Many CIT algorithms have been proposed to overcome the unstable ill-posed problem, which can be generally categorized into two groups: non-iterative and iterative algorithms.
Non-iterative tomographic algorithms for imaging the ionosphere include the regularization method [26], the singular value decomposition (SVD) method [27][28][29], the orthogonal function method [30], Bayesian statistical ionospheric tomography [31,32], ionosphere tomography based on compressed sensing [33], and so on. In recent years, there have been some impressive reports of the imaging of three-dimensional ionospheric electron density [31][32][33][34][35][36]. A novel tomography method based on Bayesian statistical inversion with prior distribution has been proposed [31,32], and Gaussian Markov random fields have been used to construct the prior electron density distribution in ionospheric 3-D multi-instrument tomography [35,36]. Furthermore, ionospheric tomography based on compressed sensing has been developed to reconstruct electron densities using a limited number of observations in certain conditions [33].
Iterative tomographic algorithms include the algebraic reconstruction technique (ART) [1] and other modified versions based on ART algorithms, such as the multiplicative ART (MART) [37,38], the simultaneous iteration reconstruction technique (SIRT) [39,40], and the simultaneous algebraic reconstruction technique (SART) [41]. Due to their simple implementation and high computational efficiency, iterative algorithms have been used in many experimental ionospheric tomography applications. However, there are also some limitations. Accurate initial values of electron density are required for these iterative reconstruction algorithms, and such values are not always straightforward to obtain. In the meantime, the limited amount of observation data used for ionospheric tomography may result in incomplete coverage of the target region. In order to improve the reconstruction quality of ionospheric electron density, a series of improved tomographic methods have been investigated. A weighted ART algorithm has been presented, with this variant of the ART algorithm assigning a weight to the density corrections that was proportional to the electron density normalized to the peak density along the GNSS rays [42,43]. Another improved ART algorithm (IART) based on the ART algorithm was proposed by introducing the electron density to each iteration procedure to compute the offset of the current electron density of the ART [44]. A multiscale tomographic model with various voxel sizes was proposed to provide a better reconstruction quality of the electron densities [45,46]. Another model with variable voxel size was also introduced [47]. Moreover, the additional information can be added using spatial-temporal constraints to improve the reliability of ionospheric tomography. A constrained SART method has been presented, which demonstrated a convergence speed significantly higher than that of the classical SART method [48]. The smoothing constraints approach of ionospheric tomography can effectively constrain near voxels with no ray path information, and a SIRT method has been proposed in which Sobolev's norm was used as a stabilizer [23]. The Phillips smoothing method and the second-order Laplace operator have been used to provide constraints regarding the smoothness of neighboring voxels [49,50]. In addition, the relaxation parameter is also essential for the ionospheric tomography, since relaxation parameter determines the degree of dependence on the initial model. Generally, the relaxation parameter has been determined on the basis of experience. In order to optimize the relaxation parameter, automatic search technology has been proposed to train the relaxation parameter using the root mean square error (RMSE) of STEC at each step of the iteration [51].
In this paper, we focus on the SART iterative reconstruction algorithms of tomography. The differences between the observed GNSS TEC and the estimated TEC are distributed among the voxels along the ray path from the satellite to the receiver in proportion to the intersection lengths of each ray-voxel of the tomographic grid. The correction assigned to a voxel of electron density is obtained by dividing the TEC difference by the proportion of the intersection length within the specific voxel in the iterative inversion process. In other words, TEC differences are redistributed among the voxels according to the lengths of the ray-voxel intersections, but the intersection length only represents the geometric contribution of the voxel to GNSS TEC. However, TEC is the integral of all electron densities along the satellite signal path, and the contribution of a single voxel to GNSS TEC includes the geometry contribution of the intersection lengths of the ray-voxel and the contribution of the electron density of that voxel. Furthermore, the inversion errors of electron densities retrieved using the GNSS tomographic reconstruction technique are the dominant source Remote Sens. 2023, 15, 2939 3 of 18 of discrepancies between the measured and calculated TEC, while the intersection length only has an amplification effect on the inversion errors of electron density. The iterative reconstruction algorithm can be improved by considering the joint contributions of the intersection length and the corresponding electron density in the specific voxel to determine the assigned value of TEC difference in a specific voxel.
In this work, an extended SART (ESART) method is proposed by introducing electron density parameters into ionospheric tomography. To account for the significant variation of the IED gradient with altitude, different height intervals of tomographic voxels along the altitude of the ionosphere are applied to the new ESART method. We intend to address the effectiveness of the ESART method in the reconstruction of IED profiles. In addition, the new method is also tested in different ionospheric conditions during geomatically quiet time as well as during geomagnetic storms in order to validate its feasibility. The data used include GNSS observations and the ionosonde data over the period of 17-19 March 2015 and 10-12 January 2019. Furthermore, an inter-voxel smoothing strategy and an iteration stop criterion are adopted, in line with Stolle's suggestion [52]. The outline of this paper is arranged as below: Section 2 gives a brief introduction of the conventional SART and ESART methods. Section 3 presents the data and the evaluation methodology of our experiment. Section 4 validates the feasibility and effectiveness of the ESART method. Finally, a conclusion is given in Section 5.

Simultaneous Algebraic Reconstruction Technique (SART)
The SART method, which was proposed by Anders Andersen and Avinash Kak in 1984 [53], is a useful tomographic technique for imaging the ionosphere when observation data are limited. Compared to the standard ART, the SART method has the advantage of being able to effectively smooth noise, enabling the generation of a high-quality reconstruction, as can be seen from Equation (1) below: where the iterative step number is k, and a given projection is i. The SART technique starts with an initial approximation x 0 of electron density; then, the corrections for every electron density voxel are computed and stored in a separate array (namely, the correction array) for the first GNSS ray. Then, we move to the next GNSS ray and update the correction array until all rays have been solved. At each iterative step, all GNSS rays in each electron density voxel are computed, and the current estimate of electron density x k−1 j is refined to a new iterative x k j by adding the correction array to the image array. The vector of electron density unknowns are updated for all GNSS rays contributing to a single voxel, and are summed before the new iteration is refined. The difference between the measured TEC (d i ) and estimated TEC ∑ N j=1 a ij x k−1 j is redistributed among N voxels and the i GNSS ray path proportionally to the length of the ray-voxel intersection a ij . The relaxation factor is represented by λ, and it is usually selected to be a constant within the interval (0.0, 1.0] to reduce the influence of noise. In our study, λ = 0.5. j(j = 1, N) is the number of the rayvoxels passed by the i GNSS ray, and i(i = 1, M) is the total number of GNSS rays. Under ideal conditions, the SART method can converge on a weighted least squares solution [54].

Extended SART (ESART) Method
From Equation (1), it can be seen that the intersection length a ij determines the allocation of the TEC difference d i − ∑ N j=1 a ij x k−1 j in the ray-voxels along the GNSS ray path in the SART method. Voxels with a longer ray-voxel intersection will obtain larger correction values for electron density than voxels with smaller ray-voxel intersections. In addition, all voxels with the same intersection length will have the same correction values, even those voxels that make different contributions to the real TEC due to differences in electron density values. Since the contribution of one voxel to the TEC along the GNSS satellite-to-receiver path is the product of the intersection length with the electron density in the corresponding voxel, the electron density is one of the main sources of inversion error. Therefore, the intersection length only has an amplification effect on the inversion error of electron density. Therefore, the differences between the TEC measured using GNSS and the TEC estimated via GNSS tomography should be redistributed along the ray-voxels while accounting for the contribution of the voxel to TEC. In this work, the electron density parameter is introduced into Equation (1), and it is expected that the proposed ESART method will outperform the SART method, with the improved inversion accuracy of electron densities. The proposed iterative algorithm can be represented using Equation (2), below: The symbols in Equation (2) have the same meanings as those in Equation (1). According to Equation (2), the electron density is considered as a parameter for the correction of individual voxels, in contrast to the SART method (as expressed in Equation (1)). The extended method ensures that those voxels that make larger contributions to TEC will receive larger corrections. This can be achieved on the basis of the ratio of the contribution of each voxel a ij x k−1 j to the integration of all voxels' contributions, ∑ N j=1 a ij x k−1 , along the GNSS line-of-sight.

Data Sources and Preprocessing Strategy
In this work, we evaluated the performance of the ESART method when imaging IED variations using GNSS data during a geomatically quiet period on 10-12 January 2019, when the Disturbed Storm Time Index (Dst) did not exceed the absolute Dst 10 nT, and an intense geomagnetic storm period on 17-19 March, when Dst reached a minimum of −234 nT.
GNSS data from 196 reference stations belonging to the Crustal Movement Observation Network of China (CMONOC) during the periods of 10-12 January 2019 and 17-19 March 2015 were used. In addition, the ionosonde data from the Beijing (40.3 • N, 116.2 • E) and Wuhan (30.5 • N, 114.6 • E) stations were used to validate the reliability of IED retrieved from the ESART method. The ionosonde data provide an independent comparison with the electron density profiles from GNSS observations by using the tomography technique. The locations of the reference stations are shown in Figure 1.
To invert the IED, the absolute slant TEC measurements are calculated every 30 s with high accuracy using dual-frequency phase-smoothed code combination observations from the ground-based GNSS stations, and the differential code biases (DCBs) are corrected with an estimated accuracy of about 0.10 ns-0.13 ns using the two-step method developed by Li et al. [55]. In order to minimize multipath errors in GNSS data, an elevation cutoff of 20 • was considered. Based on the TEC data mentioned above, the tomographic method is applied to the inversion of IED. Additionally, the International Reference Ionosphere (IRI-2016) was selected as the background to compensate for the incomplete information in the GNSS measurements for tomography. To invert the IED, the absolute slant TEC measurements are calculated with high accuracy using dual-frequency phase-smoothed code combination ob from the ground-based GNSS stations, and the differential code biases (DCB rected with an estimated accuracy of about 0.10 ns-0.13 ns using the two-st developed by Li et al. [55]. In order to minimize multipath errors in GNSS data tion cutoff of 20° was considered. Based on the TEC data mentioned above, graphic method is applied to the inversion of IED. Additionally, the Internatio ence Ionosphere (IRI-2016) was selected as the background to compensate for plete information in the GNSS measurements for tomography. Figure    To invert the IED, the absolute slant TEC measurements ar with high accuracy using dual-frequency phase-smoothed code co from the ground-based GNSS stations, and the differential code rected with an estimated accuracy of about 0.10 ns-0.13 ns usin developed by Li et al. [55]. In order to minimize multipath errors tion cutoff of 20° was considered. Based on the TEC data ment graphic method is applied to the inversion of IED. Additionally, ence Ionosphere (IRI-2016) was selected as the background to com plete information in the GNSS measurements for tomography.

Outline of the Experiment
In this paper, the feasibility and effectiveness of the new ESART method is evaluated by comparing the IED obtained using the conventional SART method and that obtained using the IRI-2016 model. The IED observations using the three approaches are compared with ionosonde data, which are considered to be a reference in our study, and these comparisons are carried out for both quiet and disturbed ionospheric conditions. The variations in the latitude-longitude maps of the IED are also compared with those in the GNSS TEC maps during geomagnetic storms. The response of the latitude-altitude map of the IED to geomagnetic storms is also discussed. Additionally, the average root square error (RMSE) and absolute error of the reconstructed IED values are calculated and used to evaluate the performance of different reconstruction methods.
As shown in Figure 1, the reconstructed ionospheric regions cover an area from 75 • E to 135 • E in longitude and from 15 • N to 55 • N in latitude, with a height range varying from 90 to 1000 km. For both the SART and ESART methods, the horizontal resolution of each prism voxel is 2.5 • × 5 • lat/lon. The vertical resolution used for both methods is also the same. However, considering the large vertical variations of the ionospehric electron density gradient in the discretized ionospheric region, the vertical intervals of electron density voxels within the height ranges 90-210 km, 210-400 km, 400-700 km and 700-1000 km are set to 30, 10, 50 and 100 km, respectively. These height intervals are designed to be smaller for regions with higher electron density gradients (especially for the peak electron density region) than for regions with smaller gradients. In addition, it is assumed that the ionosphere tends to be quasi-stationary during the effective TEC measurement time of 10 min adopted for imaging the 3-D electron densities using GNSS data in this study [42,56]. The electron density is considered to be uniformly distributed within each voxel, and to remain unchanged within the short observation times.
In addition, an intervoxel smoothing strategy is adopted to fill the gaps in the GNSS data for those voxels not intersected by any GNSS rays and which depend on the initializing values of the IRI-2016 model. During the smoothing, the electron density of each voxel is redefined as a weighted sum of itself and the values of several neighboring voxels. Furthermore, a three-dimensional distance-weighted Gaussian-like boxcar average is also used, while the information of voxels with GNSS rays may be distributed to regions with data gaps. An iteration stop criterion is used, according to the mean reconstruction residual and root mean square of the TEC [52]. According to Figures 3 and 4, it can be clearly seen that the constructed peak electron densities (NmF2) obtained using the ESART method are generally closer to those obtained using the ionosonde than those obtained using the SART method or the IRI-2016 model. In addition, there are large discrepancies in NmF2 between the IRI-2016 model and the  Figures 3 and 4, it can be clearly seen that the constructed peak electron densities (NmF2) obtained using the ESART method are generally closer to those obtained using the ionosonde than those obtained using the SART method or the IRI-2016 model. In addition, there are large discrepancies in NmF2 between the IRI-2016 model and the ionosonde results. In Figure 3(a2,a3,c1-c3) and in Figure 4 at 4:00 UT, 7:00 UT and 10:00 UT, the NmF2 values derived using the ESART method are evidently closer to those obtaind using the ionosonde than those obtained using the SART method. This indicates that the ESART method outperforms the SART method under these conditions. In addition, it can also be seen that there are small discrepancies in the NmF2 values obtained using the two tomographic methods, while the electron density profiles obtained using the IRI-2016 model are in better agreement with the ionosonde measurements. Figures 5 and 6 show that the electron density profiles obtained using the ESART method are in better agreement with the ionosonde measurements than those of the SART method or the IRI-2016 model. According to Figure 2 Similarly, it can be seen that the NmF2 values reconstructed using the ESART method are in better agreement with the profiles obtained from the ionosondes than those obtained using the SART method or IRI-2016.

According to
These results suggest that ESART can generate more accurate IED profiles than the conventional SART method or the IRI-2016 model in both geomagnetically quiet conditions and during geomagnetic storms. This suggests that our approach, in which the tomographic electron density corrections are estimated based on the product of the intercept and electron density within voxels makes the assignment of corrections at different heights more effective. Therefore, a comparison of the results further proves that the ESART method provides better electron density reconstruction than the SART method or the IRI-2016 model.
It should be noted from Figure 5(a2,b2) and Figure 6(a2-d2) that the peak height of the F2 layer (hmF2) provided by the ESART and SART methods, as well as by the IRI-2016 model, deviate significantly from the electron density profiles obtained using the ionospheric ionosonde. This indicates that none of the three methods is able to accurately describe the fluctuation in hmF2 triggered by intense geomagnetic storms.
It should also be pointed out that there are obvious discrepancies in the topside part between the electron density profiles obtained from the ionosonde measurement and those derived using other methods. This is because the ionogram provides information for directly calculating the vertical electron density profile up to the peak of the ionospheric F2 layer [59], but does not provide direct information on the topside ionosphere. The profile above the peak is estimated using an alpha-Chapman function with a scale height derived from the profile shape around the F2 peak [60]. Therefore, in this paper, the focus is on the bottom side profile of the ionosphere up to the peak of the F2 layer.
Furthermore, the errors, as represented by RMSE, and the absolute differences of the three approaches compared to the ionosonde data are calculated and compared. The RMSE values and mean absolute error of the reconstructed IED profiles can be calculated by the following equations: where Ne ionosonde i denotes the ionospheric electron density obtained from the ionosonde measurement, Ne recon i represents the reconstructed bottom electron density derived from GNSS tomography, and n is the total number of the tomographic voxels along the ionospheric electron density profile.
Tables 1-4 list the RMSE and mean absolute error of the bottom profiles of IED at the Beijing and Wuhan stations during the geomatically quiet time and during the intense geomagnetic storm, respectively. tion, the RMSE values of ESART, SART and IRI-2016 are 0.230, 0.283, 0.264 × 10 11 el/m 3 , respectively, and the absolute errors for the three methods are 0.173, 0.233, 0.185 × 10 11 el/m 3 , respectively. For Wuhan Station, the average RMSE of ESART, SART and IRI-2016 are 0.224, 0.387, 0.427 × 10 11 el/m 3 , respectively, and the average absolute errors are 0.176, 0.329, 0.329 × 10 11 el/m 3 , respectively, for the three methods. According to the results, the electron density obtained using the ESART method is in relatively good agreement with the ionosonde data, and Table 2 shows that the electron density of the SART method is slightly better than that obtained using the IRI-2016 model. Tables 3 and 4  These results also demonstrate that the RMSE and the absolute errors of the IED reconstructed using the ESART method are smaller than those when using the SART method.
On the whole, the performance of the ESART method is better than that of the IRI-2016 model, and the ESART method demonstrates an improvement in the reconstruction of IED values compared to the SART method.
As is well known, variations in NmF2 and TEC have a relatively strong correlation during the daytime. In order to further validate the effectiveness of the ESART method, in the following section, we will discuss the features among the reconstructed IED maps obtained using the IRI-2016 model, the SART and ESART methods, and the GNSS TEC maps during the geomagnetic storm.
The TEC variation maps, which characterize the evolution of the ionospheric storm in the regions of our study, were calculated using GNSS data on 17-19 March 2015. In order to analyze the responses of the electron densities obtained using the three different methods to the negative phase storm, the temporal variations in the latitude-longitude maps of the IED at an altitude of 250 km (near the altitude of peak electron density) were reconstructed using the ESART, SART, and IRI-2016 models, respectively. Figure 7 illustrates the evolution of IED and the corresponding TEC from 01:00 UT to 10:00 UT on 18 March 2015 during the main phase of the magnetic storm. The ordinate and the abscissa are the geographic latitude and longitude. The left three columns show the IED maps obtained using IRI-2016, SART, and ESART. The color bar represents TEC in TECu. The subfigures in the fourth column are two-dimensional pseudo-color maps of GNSS TEC, and the color bar represents TEC in TECu.
In Figure 7, the storm-related ionospheric disturbances can be clearly seen, and the distinctive features of a negative ionospheric storm can be observed at 1:00 UT-10:00 UT [57,58]. From the GNSS TEC maps, it can be observed that the GNSS TECs undergo a decrease, and there is also a significant decline at 1:00-10:00 UT (9:00-18:00 LT). The GNSS TEC maps clearly depict a negative phase of the strong geomagnetic storm on 18 March 2015. However, the spatio-temporal evolution characteristics of the IED background values can be characterized using the IRI-2016 model, and no ionospheric disturbance phenomenon can be observed in the left column subfigures of Figure 7. The IED values obtained using the ESART and SART methods were obviously smaller than the background values obtained using the IRI-2016 model. In addition, the IED maps obtained using the SART and ESART methods indicated an anomalous decrease in ionospheric density, with a sharp drop in TEC values. In Figure 7, the storm-related ionospheric disturbances can be clearly seen, and the distinctive features of a negative ionospheric storm can be observed at 1:00 UT-10:00 UT [57,58]. From the GNSS TEC maps, it can be observed that the GNSS TECs undergo a decrease, and there is also a significant decline at 1:00-10:00 UT (9:00-18:00 LT). The GNSS TEC maps clearly depict a negative phase of the strong geomagnetic storm on 18 March 2015. However, the spatio-temporal evolution characteristics of the IED background values can be characterized using the IRI-2016 model, and no ionospheric disturbance phenomenon can be observed in the left column subfigures of Figure 7. The IED values obtained using the ESART and SART methods were obviously smaller than the background values obtained using the IRI-2016 model. In addition, the IED maps obtained using the SART and ESART methods indicated an anomalous decrease in ionospheric density, with a sharp drop in TEC values.
It should also be noted that the electron density depletion region of the IED maps obtained using the ESART method is more obvious compared to the IED maps obtained using the SART method. It can be seen that the variations in the IED maps obtained using the ESART method and the TEC maps were almost synchronous, which demonstrates the effectiveness of the ESART method proposed in this paper.
According to related reports [57,58], a weakly positive TEC disturbance occurred at 08:00 UT in 45 N°-60 N° zone on 17 March 2015. One hour later, the positive TEC disturbances extended to 30° N. In addition, the distinct mid-latitude zones of increased TEC were then observed at 10:00-13:00 UT. It should also be noted that the electron density depletion region of the IED maps obtained using the ESART method is more obvious compared to the IED maps obtained using the SART method. It can be seen that the variations in the IED maps obtained using the ESART method and the TEC maps were almost synchronous, which demonstrates the effectiveness of the ESART method proposed in this paper.
According to related reports [57,58], a weakly positive TEC disturbance occurred at 08:00 UT in 45 • N-60 • N zone on 17 March 2015. One hour later, the positive TEC disturbances extended to 30 • N. In addition, the distinct mid-latitude zones of increased TEC were then observed at 10:00-13:00 UT.  [57,58], demonstrating the ability of the proposed ESART method to perform imaging of the ionosphere.
It should be noted that fluctuations in the hmF2 cannot be observed in Figure 8, which is consistent with the phenomenon reflected in Figure 5(a2,b2) and Figure 6(a2-d2). In other words, the ESART method is not able to reveal the fluctuations in hmF2 triggered by intense geomagnetic storms, which is also an issue that needs to be further studied in future. the reports in [57] and [58], demonstrating the ability of the proposed ESART method to perform imaging of the ionosphere.
It should be noted that fluctuations in the hmF2 cannot be observed in Figure 8, which is consistent with the phenomenon reflected in Figures 5(a2,b2) and 6(a2-d2). In other words, the ESART method is not able to reveal the fluctuations in hmF2 triggered by intense geomagnetic storms, which is also an issue that needs to be further studied in future.

Conclusions
This paper proposed the ESART method for imaging three-dimensional electron densities using an ionospheric tomography technique. In the ESART method, the product of the intersection length with the electron density is introduced into the corresponding voxel, instead of the length of the ray-voxel intersection, while the differences between the measured and calculated TEC are also reasonably redistributed among the ray-voxels. This new method ensures that the actual correction of the electron density in a voxel matches the contribution of the corresponding voxel to the GNSS TEC by accounting for the contribution of a voxel to the GNSS TEC in terms of both its geometric contribution and the contribution of its electron density.
The ESART method was evaluated by comparing it to the SART method and the IRI-2016 model under geomagnetically quiet conditions on 10-11 January 2019, as well as during a geomagnetic storm period on 17-19 March 2015. Data from an ionosonde were used as a reference in this study. The results show that the constructed electron density obtained using the ESART method was generally closer to the reference data from the ionosonde than that of the SART method or the IRI-2016 model (especially for the NmF2), confirming the effectiveness of the proposed ESART method.
Furthermore, we presented latitude-longitude pictures of the IED during the intense geomagnetic storm, displaying the spatial and temporal variations, as well as the variations in latitude profile with altitude, on 17-19 March 2015, demonstrating the great potential of the proposed ESART method to provide high-quality IED reconstruction.
It should be noted that the results presented in this study, and which demonstrate the feasibility of the new method, are preliminary. The top plasma density obtained using the ESART method should be carefully evaluated using COSMIC and Swarm data, and the data span was relatively short compared to the 11-year periodic solar variation activity. Therefore, the next step of research will focus on using more data to explore the applicability of the ESART method in the GNSS application field. In addition, the theoretical background for the ESART method has not yet been properly investigated; nevertheless, the practical results suggest that the method works as expected.