Assimilating GNSS TEC with an LETKF over Yunnan, China

Abstract

A robust ionospheric model is indispensable for providing the atmospheric delay corrections for global navigation satellite system (GNSS) navigation and positioning and forecasting the space environment. The accuracy of ionospheric models is limited due to the simplified model structures. Complicated spatiotemporal variations in total electron content (TEC) biases between GNSS and international reference ionosphere (IRI) suggest a robust strategy to optimally combine GNSS and IRI TEC for high-precision modeling. In this paper, we propose a novel ionospheric data assimilation method, which is a local ensemble transform Kalman filter (LETKF), to construct an ionospheric model over Yunnan in southwestern China. We used the LETKF method to assimilate the ionospheric TEC extracted from GNSS observations in Yunnan into the IRI-2016 model. The experimental results indicate that the ionospheric data assimilation has a more pronounced improvement effect on the IRI empirical model during periods of geomagnetic quiet than during periods of geomagnetic disturbance. On quiet magnetic days, the skill score (SKS) of the assimilation is 0.60 and the root mean square error (RMSE) values before and after assimilation are 5.08 TECU and 2.02 TECU, respectively. The correlation coefficient after assimilation increases from 0.94 to 0.99. On magnetic storm days, the SKS of the assimilation is 0.42 and the RMSE values before and after assimilation are 5.99 TECU and 3.46 TECU, respectively. The correlation coefficient after assimilation increases from 0.98 to 0.99. The results suggest that the LETKF algorithm can be considered an effective method for ionospheric data assimilation.

1. Introduction

In the near-Earth space environment, the ionosphere has an important impact on human activities as a medium for radio wave propagation. Dramatic spatial and temporal variations in the ionosphere can seriously reduce the quality of global short-wave communications and the accuracy of satellite navigation and positioning [1,2]. The ionospheric total electron content (TEC) is one of the important parameters used to characterize the spatial and temporal variability of the ionosphere [3]. The question of how to accurately obtain the characteristics and variation patterns of ionospheric TEC distribution and establish a system that both contains intrinsic physics and reflects real observations while satisfying the real-time monitoring of the space environment and providing high-precision quasi-real-time ionospheric error correction information has become one of the key focuses of research in the study of the ionosphere.

In recent years, ionospheric data assimilation has become an important method for obtaining high-precision ionospheric TEC with the improvement of the global navigation satellite system (GNSS) and the increase in the number of ground-based observation networks [4,5,6,7,8,9,10,11,12,13]. The commonly used assimilation methods are mainly divided into two categories: the variational method and the Kalman filter method. For the variational method, Bust et al. [14] used a three-dimensional variational assimilation technique to assimilate GNSS TEC into an empirical ionospheric model and constructed an ionospheric data assimilation three-dimensional (lDA3D) model, which can handle a wide range of observations such as of the altimeter, ground-based GNSS, and occultation data. The model can process a variety of observations such as of the altimeter, ground-based GNSS, and occultation data to obtain information on ionospheric variability in a uniform global grid. Mengist et al. [15,16] investigated the performance of the four-dimensional technique for ionospheric data assimilation using the empirical model of the international reference ionosphere (IRI) as a background model in order to develop a regional ionospheric forecast model. Performance was improved in regions with no observational data and also improved by adding multi-source data. The system could operate in both calm and storm conditions. Jeong et al. [17] evaluated the performance of the data assimilation model (IDA4D) by introducing multiple types of data into the empirical model around the Korean Peninsula and found that the IDA4D model can be used as a reliable nowcast model. For the Kalman filter method, Schunk et al. [18] developed a data assimilation model, the Global Assimilation of Ionospheric Measurements (GAIM), based on a physical model of the ionospheric plasma, that uses Kalman filtering as the data assimilation algorithm. Lin et al. [19] proposed a Gauss–Markov Kalman filter method to assimilate GPS receiver and space radio occultation instrument observations into the IRI model. It was shown that this method could better improve the accuracy of the data assimilation analysis, and it was confirmed that this method could be used to reconstruct the three-dimensional ionospheric electron density. Yue et al. [20] used the ensemble Kalman filter (EnKF) algorithm to assimilate the electron density data observed by the Millstone Hill incoherent scattering radar. The results showed that the EnKF technique outperformed the 3DVAR technique, especially in the region of data gaps. He et al. [21] used the EnKF algorithm to evaluate the effectivenesses of various types of information for the nowcasting and forecasting of ionospheric parameters in China and adjacent areas and found that the algorithm can be considered a useful tool for the accurate specification of ionospheric nowcasts and predictions. Kosary et al. [22] presented a sequential calibration approach, based on the EnKF needed to calibrate model parameters, that uses a short period of GNSS network measurements to improve TEC estimations. Tang et al. [23] presented a data assimilation model for a regional ionosphere based on the local EnKF and found that the method could allow effective regional ionospheric data assimilation. Forootan et al. [24] presented a data assimilation method based on principal component analysis and the EnKF. The method can be used to improve VTEC estimates for globally available ionospheric models using IGS GIM products and can be easily applied to regional data assimilation by varying the background model and the domain of the observed field. However, the deviation between IRI and GNSS is not very clear at present, and the ionosphere varies greatly across different times and spaces. It is difficult to accurately represent prior noise with simple functional models.

The LETKF is an integrated square root filter that uses a low-rank estimate of the forecast covariance matrix to compute the relevant analysis [25]. The analysis is computed independently, grid point by grid point, by simultaneously absorbing nearby observations, which allows the LETKF to be implemented efficiently on a cluster of parallel computers [26]. The local ensemble Kalman filter, an early version of the LETKF, has been applied to global forecast models and has proven to be an effective data assimilation scheme in the troposphere, especially in regions with sparse observation densities [27]. Yunnan is located in a low latitude region, where the ionosphere overhead is very active. Establishing a high-precision ionospheric assimilation model is beneficial to improve the performance of satellite navigation and positioning. In this paper, we present the LETKF assimilation method to assimilate GNSS ionospheric TEC observations into the IRI-2016 model over Yunnan in China. Observations from 27 ground-based GNSS stations are applied to construct the Yunnan regional ionospheric TEC assimilation model. Using the self-modelled data, a preliminary assessment of the assimilation results of the ionospheric TEC data from day 121 to day 151 in 2017 is carried out.

2. Local Ensemble Transform Kalman Filter Assimilation Algorithm

Since the prediction error covariance matrix is not explicitly computed in the ETKF method, but is implemented by passing the prediction ensemble perturbation matrix, the matrix dimension will be inconsistent when using the Covariance Localization (CL) method in the ETKF and performing the Schur product operation [28]. To solve this problem, the CL method can approximate the Schur product operation between the prediction error covariance matrix and the localization function by requiring a Schur product operation in the square root matrix of the prediction ensemble perturbation matrix and the localization function in the ETKF [29]. In this way, the CL method can also be used in the ETKF method for localization.

The LETKF, like the LEnKF method, also roughly consists of three processes: prediction, updating, and localization. Since the background model used in this experiment was the empirical model IRI model, as in Experiments 1 and 2, the Gauss–Markov method was used for temporal updating, and the predicted value at the current moment $t$ was calculated from the background and analytical values at moment $t-1$.

The LETKF prediction process is as follows:

$$X_t^f=X_t^b+\left(X_{t-1}^a-X_{t-1}^b\right)e^{-\Delta t/\tau }=\left(x_1^f,x_2^f,\cdots ,x_N^f\right)$$

(1)

$${\overline{x}}^f=\frac{1}{N}{\displaystyle \sum _{i=1}^N{x}_i^f}$$

(2)

$$X^{\prime f}=\frac{1}{\sqrt{N-1}}\left(x_1^f-{\overline{x}}^f,x_2^f-{\overline{x}}^f,\cdots ,x_N^f-{\overline{x}}^f\right)=\left(u_1,u_2,\cdots ,u_N\right)$$

(3)

$$P^f=X^{\prime f}{\left(X^{\prime f}\right)}^T$$

(4)

where $X_t^f$ is the set of predictions at the current moment. $X_{}^b$ is given by the background model IRI. $X_{t-1}^a$ is the set of state analysis at the previous moment. $N$ denotes the number of ensemble members. ${\overline{x}}^f$ denotes the average state of $x_i^f$ ($i=1\sim N$). $X^{\prime f}$ is the prediction ensemble perturbation matrix. $P^f$ is the prediction error covariance matrix.

The LETKF localization process is as follows:

$${\rho }_{ij}=e^{-\alpha \left({\phi }_{ij}^2/L_x^2\right)}{e}^{-\alpha \left({\theta }_{ij}^2/L_y^2\right)}$$

(5)

$$\rho =W{W}^{\mathrm{T}}=\left(w_1,w_2,\cdots ,w_L\right){\left(w_1,w_2,\cdots ,w_L\right)}^T$$

(6)

$$W={\rho }_{eigenvectors}{\left({\rho }_{eigenvalue}\right)}^{1/2}$$

(7)

$$P_{loc}^f=P^f°\rho =Z^f{\left(Z^f\right)}^{\mathrm{T}}$$

(8)

$$Z^f=W{X}^{\prime f}=\left[\left(w_1\circ u_1,w_1\circ u_2,\cdots ,w_1\circ u_N\right),\cdots ,\left(w_L\circ u_1,w_L\circ u_2,\cdots ,w_L\circ u_N\right)\right]$$

(9)

$$v_k^f={\overline{x}}^f+\sqrt{M}{z}_k,\hspace{1em}k=1,2,\cdots ,M,\hspace{1em}M=N\times L$$

(10)

$$V^f=\left(v_1^f,v_2^f,\cdots ,v_M^f\right)$$

(11)

$$P_{loc}^f=\frac{1}{M}{\displaystyle \sum _{i=1}^M\left(v_i^f-{\overline{v}}^f\right)}{\left(v_i^f-{\overline{v}}^f\right)}^{\mathrm{T}}=Z^f{\left(Z^f\right)}^{\mathrm{T}}$$

(12)

Here, ${\rho }_{ij}$ is the correlation coefficient between point $i$ and point $j$ and $\alpha$ is the distance correlation coefficient. ${\phi }_{ij}$ and ${\theta }_{ij}$ are the distances between two points in the longitude and dimensional directions, respectively. $L_x$ and $L_y$ denote the correlation distances in the longitude and dimensional directions, respectively. In this paper, the correlation coefficient $e^{-\alpha }$ is taken to be 0.75, the correlation distance in the dimensional direction is $10°$, and the correlation distance in the longitude direction varies from about $40°$ in mid-latitudes to about $20°$ in the equatorial region. $W$ denotes the square root matrix of $\rho$, which is obtained by decomposing the eigenvalues of $\rho$ and arranging the obtained eigenvalues in descending order, taking the first 10 eigenvalues and the eigenvector pair composition; $L$ is taken as 10. $P_{loc}^f$ is the prediction error covariance matrix after localization by the CL method. $Z^f$ denotes the square root matrix of $P_{loc}^f$, obtained by using the ensemble expansion technique, with a matrix size of $n\times \left(N\times L\right)$ where $n$ denotes the number of state values in a single ensemble. $z_k$ denotes the $k$ column of matrix $Z^f$. $M$ is the size of the set after localization expansion. $V^f$ is the reconstructed prediction ensemble after localization. Equation (11) gives the reconstructed prediction error covariance matrix $P_{loc}^f$ and its square root matrix $Z^f$.

The LETKF update process is as follows:

$${\widehat{Y}}^f=R^{-\frac{1}{2}}H{Z}^f$$

(13)

$${\widehat{Y}}^f=U\Sigma V^T$$

(14)

$$Z^a=Z^{f}U{\left(\Sigma {\Sigma }^T+I\right)}^{-\frac{1}{2}}{U}^T$$

(15)

$${\overline{v}}^a={\overline{v}}^f+Z^{f}U\Sigma {\left({\Sigma }^T\Sigma +I\right)}^{-1}{V}^T{R}^{-\frac{1}{2}}\left[y-H\left({\overline{v}}^f\right)\right]$$

(16)

$$X^a={\overline{v}}^a{1}_N^T+\sqrt{N-1}{Z}^{a}randn\left(M,N\right)$$

(17)

Here, ${\widehat{Y}}^f$ is the annotated prediction–observation ensemble perturbation matrix. $R$ is the observation error covariance matrix. $U$ and $V$ are orthogonal matrices and $\Sigma$ is the singular value diagonal matrix of ${\widehat{Y}}^f$; all three are derived from a singular value decomposition of ${\widehat{Y}}^f$. $Z^a$ is the square root matrix of the analysis error covariance matrix. ${\overline{v}}^a$ is the mean of the assimilation analysis set. $X^a$ is the assimilation analysis ensemble. $1_N^T$ denotes an $1\times N$ dimensional row vector in which all elements are 1s. $randn\left(M,N\right)$ denotes an $M\times N$ dimensional random matrix with elements obeying the standard normal distribution, which completes the selection of the $N$ column analysis set; the matrix $Z^a$ is randomly transformed into an $n\times N$ dimensional matrix. The implementation steps are based on those used by Hunt et al. [25].

3. Experimental Data Accuracy Assessment

3.1. Experimental Data

In this experiment, data from 27 GNSS stations in the Crustal Movement Observation Network of China (CMONOC) were used for a 31-day ionospheric assimilation experiment that included the annual cumulative days from day 121 to day 151 in 2017 (i.e., 1 to 31 May 2017). The longitude of the region is from $97°\mathrm{E}$ to $107°\mathrm{E}$ and the latitude is from $21°\mathrm{N}$ to $29°\mathrm{N}$, with a latitude and longitude resolution of $0.5°\times 0.5°$. As shown in Figure 1, the black triangles indicate the GNSS station distribution.

Figure 2 shows the trends of the Dst index and Kp index during the period comprising days 121 to 151 in 2017. It can be seen that a magnetic storm occurred from day 147 to day 148 in 2017. The Dst index reached a minimum value of −125 nT at 08:00 UT on day 148 in 2017, and the Kp index reached a maximum value of 7 between 04:00 UT and 06:00 UT on day 148 in 2017.

In this experiment, it was difficult to include comparative data in the experimental analysis because of the large spatial resolution of Center for Orbit Determination in Europe (CODE) TEC data and the small distribution of Massachusetts Institute of Technology (MIT) TEC data in the region of Yunnan. In this study, TEC maps established by using the low-degree spherical harmonic function were used as true values for comparative analysis [30]. In the following sections, we refer to the TEC maps as self-modeling TECs.

3.2. Accuracy Assessment

To assess the accuracy of assimilation results, the RMSE can be used as a criterion for accuracy assessment; it is defined as follows:

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{n=1}^N{\left({\mathrm{TEC}}_n^{\mathrm{ass}}-{\mathrm{TEC}}_n^{\mathrm{rea}}\right)}^2}}$$

(18)

where $N$ is the total number, ${\mathrm{TEC}}_n^{\mathrm{ass}}$ is the TEC value at point $n$ obtained by data assimilation, and ${\mathrm{TEC}}_n^{\mathrm{rea}}$ is the truth data.

SKS is a common method for the quantitative assessment of assimilation effects, allowing a comprehensive assessment of the performance of ionospheric data assimilation. Three types of data are required for the calculation of SKS, namely the present or forecast value to be compared, the reference value, and the true value, which are calculated as follows [31]:

$${\mathrm{RMSE}}_{ass}=\mathrm{RMSE}\left({\mathrm{TEC}}_{ass}-{\mathrm{TEC}}_{truth}\right)$$

(19)

$${\mathrm{RMSE}}_{ref}=\mathrm{RMSE}\left({\mathrm{TEC}}_{ref}-{\mathrm{TEC}}_{truth}\right)$$

(20)

$$\mathrm{SKS}=1-{\mathrm{RMSE}}_{ass}{/\mathrm{RMSE}}_{ref}$$

(21)

where ${\mathrm{RMSE}}_{ass}$ and ${\mathrm{RMSE}}_{ref}$ denote the root mean square errors of the present or forecast and reference values to be compared, respectively, and ${\mathrm{TEC}}_{ass}$, ${\mathrm{TEC}}_{ref}$, and ${\mathrm{TEC}}_{truth}$ denote the present or ionospheric TEC data to be compared with the forecast, reference, and true values, respectively. From the equation, it can be derived that $\mathrm{SKS}\le 1$, and the larger $\mathrm{SKS}$ is, the better the value to be compared is. When $\mathrm{SKS}=1$, it means that the value to be compared is exactly the same as the true value. When $\mathrm{SKS}=0-1$, it indicates that the value to be compared is better than the reference value. When $\mathrm{SKS}<0$, it indicates that the value to be compared is worse than the reference value.

4. Results and Analysis

4.1. Ionospheric Geomagnetic Calm Conditions

To verify the reliability of the assimilation method, a comparison of the true TEC and the TECs before and after assimilation for two GNSS stations, KMIN and XIAG, for the 10-day period from day 125 to day 134 in 2017 is given in Figure 3. From Figure 3, it can be seen that the assimilated TEC data is closer to the true value with good agreement compared to the background model IRI data.

The daily average RMSE of the ionospheric TECs before and after the assimilation of the data and the daily average SKS of the assimilated data from day 125 to day 134 in 2017—i.e., for a 10-day period—are given in Figure 4, with the self-modeled data included as true values for comparison. As can be seen in Figure 4, the background data IRI TEC RMSE before assimilation deviates more from the true value. After assimilation, it is closer to the true value and the daily average RMSE decreases significantly. In Figure 4, all SKS values are greater than 0 during the aforementioned 10 days. This result indicates a certain improvement of IRI TEC relative to the background data.

The SKS values of the data-assimilated ionospheric TEC and the RMSEs before and after assimilation over time are given in Figure 5 for 24 h on day 126, 2017. It can be seen from Figure 5 that after data assimilation, the SKS values of the assimilated TEC are greater than 0, which indicates that after data assimilation, the accuracy of the assimilated TEC has a significant improvement compared to the IRI TEC. The greatest improvement relative to the IRI model for this day occurred at 06:00 UT with an SKS of 0.92 and RMSEs of 8.19 TECU and 0.67 TECU for IRI TEC and the assimilated TEC, respectively. The average SKS for this day was 0.72 and the average RMSE for the TECs before and after assimilation was 4.19 TECU and 1.18 TECU, respectively.

Figure 6 shows the distribution of ionospheric TEC data in Yunnan and the surrounding regions for ionospheric pierce point (IPP) TEC, background model IRI TEC, assimilated TEC, and true TEC for the observed data at 09:00 UT on day 132, 2017. As can be seen in Figure 6, overall, the data in the background model IRI TEC is significantly overestimated compared to the true value while the assimilated TEC is closer to the IRI TEC than the IRI TEC is to both the true and IPP TECs. This suggests that the assimilated TEC has better agreement with the true value compared to the IRI TEC.

Figure 7 gives the TECs before and after data assimilation compared to the true data on day 127 in 2017. The histograms of the residuals of the background model IRI TEC and assimilated TEC compared to the true values are given in Figure 7a,b, respectively, which show that the residual distribution of the assimilated TEC data is closer to the unbiased Gaussian distribution compared to the IRI TEC. The mean and RMSE values of the residuals between the IRI TEC and the true value are −2.16 TECU and 5.08 TECU, respectively, while the mean and RMSE values of the residuals between the assimilated TEC and the true value drop to 0.25 TECU and 2.02 TECU, respectively. Figure 7c,d give scatter plots of the TECs before and after assimilation compared with the true value, respectively, and it can be seen that the correlation coefficient with the true value increases from 0.94 to 0.99 before assimilation.

Figure 8 shows the background model IRI TEC, true TEC, and assimilated TEC distributions for ionospheric TEC in Yunnan and the surrounding regions for 132 full days in 2017. As can be seen in Figure 8, the background model gives an obvious overestimation of the ionospheric TEC compared to the true TEC distribution, especially at the 06:00 UT and 10:00 UT moments, where the difference is obvious, with a maximum difference of 12 TECU, while the TEC distribution after assimilation is more consistent with the true TEC distribution.

4.2. Under Ionospheric Geomagnetic Disturbance Conditions

Figure 9 gives a comparison of the true TEC and the TECs before and after assimilation for the KMIN and XIAG observatories during the 10-day period from day 142 to 151 in 2017; a geomagnetic disturbance occurred on day 148. As can be seen in Figure 9, similar to in Figure 3, the agreement between assimilated TEC and the true value is significantly improved compared to IRI TEC on geomagnetically calm days while the agreement between assimilated TEC and the true value is weakened on the perturbed days 147 to 148 but still improved compared to IRI TEC.

The daily average SKS of the data assimilation and the daily average RMSE of the ionospheric TECs before and after the assimilation are given in Figure 10 for days 142 to 148 during the 10-day period in 2017. It can be seen in Figure 10 that, similar to in Figure 4, the RMSE values of the assimilated TEC during the 10-day period are reduced compared to the pre-assimilation IRI TEC, and the SKS values are all greater than 0, indicating that the assimilated TEC data are improved relative to the background model IRI TEC data. It can also be seen that the SKS of the assimilated TEC on the 148th geomagnetic disturbance day is reduced compared to the calm day, i.e., the difference between the RMSE of the assimilated TEC and the IRI TEC is reduced.

The SKS of the ionospheric TEC assimilated by the data and the RMSE values before and after assimilation over the period of 24 h during the geomagnetic disturbance event day in 2017, day 148, are given in Figure 11. As can be seen in Figure 11, the assimilated SKS is greater than 0 for most moments, and only a few moments show negative values, i.e., the RMSE of the assimilated TEC is greater than that of the IRI TEC, and the assimilation is not as effective as the background data before assimilation. The maximum and minimum values of SKS during this day were 0.85 and −0.21, respectively. The average SKS was 0.42, and the average RMSE values before and after assimilation were 5.99 and 3.46, respectively. From the daily average values of SKS and RMSE, we conclude that the data assimilation still had some improvement effect on the background model on the geomagnetic disturbance day.

Figure 12 shows the distribution of IPP TEC, the background model IRI TEC, assimilated TEC, and true TEC for the observed data at 06:00 UT on day 148 in 2017, the geomagnetic disturbance day, for the ionospheric TEC distribution in Yunnan and the surrounding region. As can be seen in Figure 12, the values of the background model IRI TEC are significantly lower compared to the true value, especially between 25°N and 29°N. The assimilated TEC is closer to the real value. This shows that the assimilated TEC is in better agreement with the true value than the IRI TEC on stormy days.

Figure 13 shows the TECs before and after data assimilation compared to the true data for day 148 in 2017, the day of geomagnetic disturbance. The residual histograms of the background model IRI TEC and assimilated TEC compared to the true value are given in Figure 13a,b, respectively, which show that the residual distribution of the assimilated TEC is closer to the unbiased Gaussian distribution compared to the IRI TEC. The mean and RMSE values of the residuals between the IRI TEC and the true value are −5.32 TECU and 5.99 TECU, respectively, while the mean and RMSE values of the residuals between the assimilated TEC and the true value drop to −3.22 TECU and 3.46 TECU, respectively. Figure 13c,d give the scatter plots of TEC before and after assimilation compared with the true value, respectively, and it can be seen that the correlation coefficient with the true value increases from 0.98 to 0.99 before assimilation.

Figure 14 shows the background model IRI TEC, true TEC, and assimilated TEC distributions of ionospheric TEC in Yunnan and the surrounding region throughout the day on the geomagnetic disturbance day in 2017, day 148. It can be seen from Figure 14 that the background model gives a larger difference between the ionospheric TEC and the true value compared to the true TEC distribution, especially at 10:00 UT, 14:00 UT, and 18:00 UT, where the IRI TEC is significantly underestimated, while the TEC distribution after assimilation is not very different from the true value but is not as effective as the assimilation distribution on geomagnetically calm days.

5. Discussion

To make for a better statistical analysis, Figure 15 shows the average monthly RMSE values of TEC before and after assimilation in 2017. The RMSE after assimilation is larger in May and September 2017. The poor performance of these two months is due to the magnetic storms that occurred in both months. In general, the assimilation method performs well on different days, months, and times of the day under different geomagnetic conditions.

The LETKF is an integrated square root filter that uses a low-rank estimate of the forecast covariance matrix to compute its analysis. This analysis is computed independently, grid point by grid point, by simultaneously absorbing nearby observations. First, the background perturbation matrix in the model and observation space is computed globally. Then, for each grid point, the observation vector mapped to the observation space and the associated background perturbation matrix are restricted to include only observations from the region surrounding the point [26]. Under quiet geomagnetic conditions, the residual distribution before data assimilation is significantly eccentric and deviates from the self-modeled TEC (Figure 7). Under disturbed geomagnetic conditions, the residual distribution before data assimilation is significantly eccentric and deviates from the self-modeled TEC (Figure 13). In Figure 4, the average daily SKS of assimilated TEC is 0.66 and the average daily RMSE values of TEC before and after assimilation during the 10 days from day 125 to day 134 in 2017 are 5.33 TECU and 1.69 TECU, respectively. In Figure 10, the average daily SKS of assimilated TEC is 0.64 and the average daily RMSE values of TEC before and after assimilation during the 10 days from day 142 to day 151 in 2017 are 4.79 TECU and 1.68 TECU, respectively. For the 31 days from day 121 to day 151 in 2017, the average daily SKS of assimilated TEC is 0.55 and the average daily RMSE values of TEC before and after assimilation are 5.03 TECU and 2.13 TECU, respectively. These experimental results show a generally good agreement between the self-modeled TEC and the assimilated TEC. The LETKF model can be effectively used for ionospheric data assimilation.

6. Conclusions

In this paper, the first ionospheric TEC assimilation modeling using the LETKF assimilation method was conducted to establish a regional ionospheric TEC assimilation model for Yunnan using the IRI model as the background model and the observation data from 27 ground-based GNSS stations. The results show that the improvement effect of data assimilation on the IRI background model is more obvious in periods of geomagnetic calm than in periods of geomagnetic disturbance. During the quiet period on day 127 in 2017, the SKS after assimilation is 0.60, the RMSE values before and after assimilation are 5.08 TECU and 2.02 TECU, respectively, and the correlation coefficients after assimilation are 0.94 and 0.99. During the storm period on day 148 in 2017, the SKS after assimilation is 0.42, the RMSE values before and after assimilation are 5.99 TECU and 3.46 TECU, respectively, and the correlation coefficients after assimilation are 0.98 and 0.99. The assimilation method performs well on different days, months, and times of the day under different geomagnetic conditions. The data assimilation can effectively fuse the input observations in both geomagnetic calm and geomagnetic disturbance periods, and the assimilated TEC accuracy is significantly improved compared with the pre-assimilation IRI TEC.