Ionospheric Assimilation of GNSS TEC into IRI Model Using a Local Ensemble Kalman Filter

: Ionospheric total electron content (TEC) is important data for ionospheric morphology, and also an important parameter for ionospheric correction in Global Navigation Satellite System (GNSS) precise positioning, navigation, and radio science. In this study, we present a data assimilation model for regional ionosphere based on a local ensemble Kalman ﬁlter (LEnKF) with the International Reference Ionosphere 2016 (IRI-2016) model as the background, to assimilate ionospheric TEC observations from GNSS. To demonstrate the results, the TEC estimates from the Crustal Movement Observation Network of China (CMONOC), which is about 260 stations in China, are applied as observation. The assessments are performed against the TEC estimates from BeiDou Navigation Satellite System (BDS) geostationary earth orbit (GEO) and against the ﬁnal products from the Center for Orbit Determination in Europe (CODE). The assimilation results are in good agreement with BDS GEO TEC and the CODE TEC on a quiet or disturbed day. The correlation coefﬁcient after assimilation is increased by about 17% compared with that before assimilation, and the RMSE after assimilation is decreased by about 42% compared with that before assimilation. Furthermore, the assimilated method is also evaluated in the single-frequency precise point positioning (PPP). The experimental results indicate that the PPP/Assimilated method can improve the GNSS positioning accuracy more effectively in comparison to the PPP/CODE. These results reveal that the LEnKF method can be considered as a useful tool for ionospheric assimilation.


Introduction
Ionospheric reflection, refraction, scattering, and absorption of radio signals are an important source of errors in velocity measurement, timing, navigation, and positioning of Global Navigation Satellite System (GNSS). Total electron content (TEC) is an important parameter for the ionosphere. An accurate ionospheric TEC model plays a significant role in improving the accuracy of GNSS precise positioning and navigation, and understanding the structure and variation of the ionosphere. To learn about ionospheric weather and improve the accuracy of GNSS navigation and positioning, numerous empirical, physics-based theoretical and mathematical models of the ionosphere have been developed. Ionospheric data assimilation is the process of combining observations with a model of the ionosphere which can effectively express the ionospheric morphology and improve the accuracy of the ionospheric TEC models to understand space weather and correct the ionospheric error in GNSS navigation and positioning. It is of great importance to provide accurate and reliable ionospheric TEC products as practically as possible [1].
In recent years, data assimilation models have been developed for global or regional ionospheric TEC models to predict ionospheric space weather [2][3][4][5][6][7][8][9][10][11]. Some powerful techniques to assimilate data into ionospheric models are variational algorithm, a Kalman filter, and an ensemble Kalman filter. Bust et al. [12] used a three-dimensional variational algorithm to assimilate GNSS TEC into an empirical ionospheric model, which incorporates available data, the associated data error covariances, a reasonable background specification, and the expected background error covariance into a coherent specification on a global grid. Aa et al. [13] developed a regional 3D ionospheric electron density model over China and adjacent areas with a three-dimensional variational method. Mengist et al. [14] applied a four-dimensional variational algorithm with a fidelity that improves with inclusion of multiple data types to a regional ionospheric nowcast system during quiet and storm days. Wang et al. [15] used a three-dimensional variational algorithm to derive global ionosphere maps (GIM) using GNSS measurements. Compared with the variational algorithm, the advantage of the Kalman filter is that it can dynamically propagate the background error covariance forward [16]. Schunk et al. [17] developed a physics-based data assimilation model-the global assimilation of ionospheric measurements (GAIM)-which would use a physics-based ionosphere-plasmasphere model and a Kalman filter as a basis for assimilating a diverse set of measurements. Aa et al.
[1] used a Kalman filter method to implement the data assimilation with GNSS data, and demonstrates the effectiveness of this method in providing accurate regional distribution of the ionospheric TEC. Lin et al. [18] presented a new ionospheric data assimilation procedure-the Gauss-Markov Kalman filter-to assimilate ground-and space-based TEC. Pasumarthi and Devanaboyina [19] used the Kalman filter method to generate hourly assimilated Indian regional TEC maps by the process of data assimilation. Qiao et al. [20] used the Gauss-Markov Kalman filter to develop a TEC model in China and adjacent regions. Compared with the Kalman filter algorithm, the advantage of the ensemble Kalman filter (EnKF) is that the error statistics are calculated from an ensemble of the forward model forecasts which run in parallel [21]. Yue et al. [21] used an EnKF technique to assimilate a one-dimensional midlatitude ionospheric theoretical model and found this technique had a better performance than the three-dimensional variational technique. Chartier et al. [22] assimilated TEC observations into a couped thermosphere-ionosphere model using an EnKF during storms, which showed the potential for the EnKF method to improve midlatitude storm-time TEC forecasting efforts. Chen et al. [23] used an EnKF method to assimilate ground-based GPS TEC observations into a theoretical numerical model of the thermosphere and ionosphere. He et al. [6] evaluated the quasi-realistic ionosphere forecasting capability by an EnKF ionosphere and thermosphere data assimilation algorithm. He et al. [24] also developed an EnKF ionosphere and thermosphere data assimilation model based on different observation systems over China and found that this model can greatly improve the quality of ionosphere specification. Although there are many good ionospheric assimilation models, it is of great significance to continue to study ionospheric assimilation methods and practical applications. At present, many data assimilation models have been developed for global or regional ionospheric TEC models to predict ionospheric space weather. However, few studies have used GNSS positioning to evaluate the accuracy of ionospheric assimilation models. Although some outstanding ionospheric data assimilation products have been developed in China, there are limited studies of data assimilation and application in China seeking to generate a regional TEC product with high precision. As an empirical model, the International Reference Ionosphere (IRI) model, which is an international project, plays an important role in the Earth's ionospheric studies and applications [25,26]. The latest version of the IRI model is IRI-2016 [27].
In this study, we propose a local ensemble Kalman filter (LEnKF) to assimilate GNSS ionospheric TEC observations into the IRI-2016 model. The TEC estimates from the Crustal Movement Observation Network of China (CMONOC), which is about 260 stations in China, are applied as observation. The assessments are performed against the TEC estimates from BeiDou Navigation Satellite System (BDS) geostationary earth orbit (GEO) and against the final products from the Center for Orbit Determination in Europe (CODE). In addition, we use the single-frequency precise point positioning (PPP) to evaluate the accuracy of ionospheric assimilation models for the first time.

Methodology
EnKF uses sample statistics to calculate the background error covariance, which is prone to negative problems such as filter divergence, underestimation of the prediction error covariance, and spurious correlation. As a result of these problems, EnKF may produce a suboptimal analysis result. In order to overcome the negative effects of undersampling, covariance localization should be adopted. The local ensemble Kalman filter (LEnKF) data assimilation method has been described in detail and mathematically derived in [28]. In this section, we only provide a brief derivation of this method to understand the implementation of the scheme on the ionospheric data assimilation.

Ensemble Kalman Filter
We all know that the ensemble Kalman filtering algorithm is developed from the Kalman filtering algorithm, and the major difference between them is that the background error covariance matrix of the ensemble Kalman filtering is statistically derived from a finite number (usually 10 to several hundred) of samples obtained from a fully nonlinear forecast. These ensemble forecasts start with a set of initial conditions, a short period of forecasting followed by an analysis process, and then continue the cycle of forward forecasting. Thus, the ensemble Kalman filter is an approximation and extension of the Kalman filter [29]. It can be written as where X a is the analysis value, X f is the ensemble TEC of the background given by the IRI model, Y obs is the TEC of GNSS observations, K is the Kalman gain matrix. P f is the background error covariance, H is the observation operator, R is the observation error covariance, X f denotes the ensemble mean TEC of the background model, Y obs denotes the ensemble mean TEC of GNSS observations, N is ensemble number. H is to make the observation field and the background have the same dimension. It is usually a matrix with different number of rows and columns at different times. In this paper, the number of rows is equal to the number of observed grid points, and the number of columns is equal to the number of grid points in the background field. The background field in the next step can be updated in time by a Gauss-Markov method, which assumes that the correction of the analyzed values to the background values decays exponentially with time and can be expressed as: where X f t+1 is the next moment forecast value, ∆t is the time step, and τ denotes the ionospheric time dependent scale. The relevant scale factor can be written as τ = 1/4 * cos(π/12 * (t − 14)) + 3/4 [20].
The ensemble-based background error covariance P f is derived from the statistics of the ensemble samples, which are initialized by randomly perturbing the F10.7 parameter input of IRI-2016 model. By calculating the standard deviation of F10.7, the error set that satisfies the Gaussian distribution mean of 0 is obtained, which is used as the initial condition to drive the IRI-2016 model.

Local Method
To eliminate the pseudo correlation, which is far from the observation, covariance localization is used, which implements the Schur product operation between the prediction error covariance matrix and the localization function to solve the spurious correlation problem by truncating the error correlation beyond a pre-specified distance in the prediction error covariance matrix. In addition, the Schur product operation can improve the rank of the prediction error covariance matrix.
It is generally believed that the ionospheric spatial correlation is Gaussian distributed, but there is no good conclusion on the distance dependence. The ionospheric distance dependence can be regarded as anisotropic. In this paper, the ionospheric distance dependence is regarded as the product of longitude, latitude and altitude directions [30]. The ionospheric correlation scales vary significantly in different directions, with significant variations in latitude, altitude, season, local time, and solar activity. The original P f is updated by the Schur product of the distance-based correlation coefficient matrix ρ and the background error covariance P f [31].
For the purpose of this paper, the correction factor can be expressed as: where ρ ij is the correlation coefficient between point i and point j, α is the distance correlation coefficient, ϕ ij and θ ij are the distances between two points in the longitude and latitude directions, respectively. L x and L y denote the correlation distances in the longitude and latitude directions, respectively. The correlation distance in the latitude direction is 10 • , and the correlation distance in the longitude direction changes linearly from about 40 • at mid-latitudes to about 20 • in the equatorial region [30]. In this paper, the correlation coefficient e −α defined at the correlation distance is 0.75, α = 0.29.

Data
The GNSS observation data from about 260 stations of the CMONOC [32] and the BDS GEO from the multi-GNSS experiment (MGEX) stations are used to investigate the ionospheric assimilation model over China and surrounding areas. Based on the observation data from CMONOC and MGEX stations, we calculate the vertical total electron content (VTEC) on all ionospheric pierce points (IPPs) by the uncombined PPP method [33], which is based on the ionosphere's single-layer model [34]. The assimilation region extends from 70 • E to 140 • E in longitude, 15 • N to 55 • N in latitude. The spatial resolution is 1 • × 1 • in the longitudinal and latitudinal directions. As shown in Figure 1, the red circles and green pentagrams denote the locations of CMONOC and MGEX stations over China and surrounding areas, respectively.
The IRI model is a joint undertaking by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI) with the goal of developing and improving an international standard for the parameters in Earth's ionosphere [35]. IRI is a widely used and recognized international reference ionospheric model. It is an empirical model based on a large number of ground and satellite observations of ionospheric parameters. The IRI model describes average conditions quite well and has shown excellent results in comparison with other models. In this paper, we use the latest version IRI-2016 as the background model. CODE GIMs are generated on a daily basis at CODE using data from about 300 GNSS sites of the International GNSS Service (IGS) and other institutions. The VTEC is modeled in a solar-geomagnetic reference frame using a spherical harmonics expansion up to degree and order 15. Piece-wise linear functions are used for representation in the time domain.
The time interval of their vertices is 1 h, which produces 1 hourly snapshots of the global ionosphere. Each TEC map has a spatial resolution of 2.5 • × 5 • in the geographic latitude and longitude.
The ionospheric data assimilation experiment is carried out by selecting the GNSS data of CMONOC and MGEX stations from 5 September 2017 to 11 September 2017. Figure  The IRI model is a joint undertaking by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI) with the goal of developing and improving an international standard for the parameters in Earth's ionosphere [35]. IRI is a widely used and recognized international reference ionospheric model. It is an empirical model based on a large number of ground and satellite observations of ionospheric parameters. The IRI model describes average conditions quite well and has shown excellent results in comparison with other models. In this paper, we use the latest version IRI-2016 as the background model. CODE GIMs are generated on a daily basis at CODE using data from about 300 GNSS sites of the International GNSS Service (IGS) and other institutions. The VTEC is modeled in a solar-geomagnetic reference frame using a spherical harmonics expansion up to degree and order 15. Piece-wise linear functions are used for representation in the time domain. The time interval of their vertices is 1 h, which produces 1 hourly snapshots of the global ionosphere. Each TEC map has a spatial resolution of 2.5° × 5° in the geographic latitude and longitude.
The ionospheric data assimilation experiment is carried out by selecting the GNSS data of CMONOC and MGEX stations from 5 September 2017 to 11 September 2017.

Accuracy Assessment
To evaluate the accuracy of the assimilation results, the root mean square error (RMSE) is used as the standard of accuracy assessment. Note that the RMSE is defined as

Accuracy Assessment
To evaluate the accuracy of the assimilation results, the root mean square error (RMSE) is used as the standard of accuracy assessment. Note that the RMSE is defined as where N is the total number, TEC rea n is the TEC value obtained by CODE or BDS GEO in the nth point, TEC rea n is the TEC value obtained by data assimilation in the nth point. CODE TEC is obtained from CODE GIM product.
We also use the kinematic single-frequency PPP to evaluate the accuracy of the ionospheric assimilation model. Three modes of kinematic PPP are performed and analyzed, which are the ionosphere-free combination model (PPP/IF), the single-frequency L1 observation based on CODE GIM constraint (PPP/CODE), and the single-frequency L1 observation based on the ionospheric assimilated model (PPP/Assimilated), respectively. The reference coordinates are obtained from the final solution by the static dual-frequency PPP [36].  Figure 3a is IPP TEC which is derived from GNSS observations. Figure 3b is IRI-2016 TEC which is derived from the IRI-2016 model. Figure 3c is assimilated TEC which is derived from the assimilated method of the local ensemble Kalman filter. Figure 3d is GIM TEC which is derived from CODE. It can be seen that the IRI-2016 TEC values are significantly overestimated in comparison to CODE TEC between 130 • E to 140 • E and 15 • N to 35 • N, and are significantly underestimated in comparison to CODE TEC between 70 • E to 140 • E and 35 • N to 55 • N. These suggest that the assimilated TEC has better agreements with CODE TEC than with the IRI-2016 TEC. To a certain extent, it indicates the improvement of the assimilation model when the GNSS-observed TEC data is assimilated into the IRI-2016 model. Figure 4 shows distributions of absolute differences and TEC correlations on 6 September 2017. Figure 4a presents absolute differences between the assimilated TEC and the CODE TEC. Figure 4b presents absolute differences between the IRI-2016 TEC and the CODE TEC. Figure 4c presents assimilated TEC distribution versus the CODE TEC. Figure 4d presents IRI-2016 TEC distribution versus the CODE TEC. Figure 4a,b present the residual histogram distributions for the data assimilation TEC and IRI-2016 TEC, respectively. It can be seen that the residual distribution before data assimilation is obviously deviated from the center and has a large deviation from CODE TEC. The residuals mean and the RMSE between the IRI-2016 TEC and the CODE TEC are −3.33 TECU and 5.20 TECU, respectively. The residual distribution after data assimilation is closer to an unbiased Gaussian distribution. The residuals mean and the RMSE between the assimilated TEC and the CODE TEC decrease to −1.60 TECU and 2.95 TECU, respectively. As can be seen from Figure 4d Figure 4a presents absolute differences between the assimilated TEC and the CODE TEC. Figure 4b presents absolute differences between the IRI-2016 TEC and the CODE TEC. Figure 4c presents assimilated TEC distribution versus the CODE TEC. Figure  4d presents IRI-2016 TEC distribution versus the CODE TEC. Figure 4a,b present the residual histogram distributions for the data assimilation TEC and IRI-2016 TEC, respectively. It can be seen that the residual distribution before data assimilation is obviously deviated from the center and has a large deviation from CODE TEC. The residuals mean and the RMSE between the IRI-2016 TEC and the CODE TEC are −3.33 TECU and 5.20 TECU, respectively. The residual distribution after data assimilation is closer to an unbiased Gaussian distribution. The residuals mean and the RMSE between the assimilated TEC and the CODE TEC decrease to −1.60 TECU and 2.95 TECU, respectively. As can be seen from Figure 4d, correlation coefficients of 0.97 and 0.95 are obtained for assimilated TEC and IRI-2016 TEC, respectively. The assimilation results are in good agreement with the CODE TEC.  To analyze the performance of TEC, we use three modes of kinematic PPP which are PPP/IF, PPP/CODE, PPP/Assimilated. Figure 5 shows the RMSE of the 3D coordinate of three modes with each epoch at different stations for day 249, 2017 (6 September 2017). The RMSE of the 3D is computed by taking the RMSE of the 3D coordinate (North, East, Up component) solutions. PPP/IF is PPP of the ionospheric-free observation which is used by dual-frequency combination method. PPP/Assimilated is PPP of the single-frequency observation based on the ionospheric assimilated TEC. PPP/IF can be used as reference accuracy. If the positioning accuracy of the single-frequency PPP with ionospheric constraints is slightly lower than or equal to that of PPP/IF, it can be proven that the iono- To analyze the performance of TEC, we use three modes of kinematic PPP which are PPP/IF, PPP/CODE, PPP/Assimilated. Figure 5 shows the RMSE of the 3D coordinate of three modes with each epoch at different stations for day 249, 2017 (6 September 2017). The RMSE of the 3D is computed by taking the RMSE of the 3D coordinate (North, East, Up component) solutions. PPP/IF is PPP of the ionospheric-free observation which is used by dual-frequency combination method. PPP/Assimilated is PPP of the singlefrequency observation based on the ionospheric assimilated TEC. PPP/IF can be used as reference accuracy. If the positioning accuracy of the single-frequency PPP with ionospheric constraints is slightly lower than or equal to that of PPP/IF, it can be proven that the ionospheric TEC is more accurate. It can be seen that the PPP/Assimilated and PPP/CODE solution present a good performance as well as PPP/IF for all stations. The results indicate that ionospheric constraint of PPP/Assimilated can effectively perform single-frequency PPP positioning. In Figure 5d, the performances of all PPP methods are generally worse due to the more active ionospheric activity and the higher values of TEC at low latitudes.   Figure 6a is IPP TEC, which is derived from GNSS observations. Figure 6b is IRI-2016 TEC, which is derived from the IRI-2016 model. Figure 6c is the assimilated TEC, which is derived from assimilated method of the local ensemble Kalman filter. Figure 6d is GIM TEC, which is derived from CODE. It can be seen that the IRI-2016 TEC values are significantly underestimated in comparison to CODE TEC between 70°E to 140°E and 15°N to 35°N. As can be seen in Figure 6b-d, the assimilated TEC has better agreements with CODE TEC than with the IRI-2016 TEC. To a certain extent, it also indicates the improvement of the assimilation model when the GNSS-observed TEC data is assimilated into the IRI-2016 model on the disturbed day.   Figure 6a is IPP TEC, which is derived from GNSS observations. Figure 6b is IRI-2016 TEC, which is derived from the IRI-2016 model. Figure 6c is the assimilated TEC, which is derived from assimilated method of the local ensemble Kalman filter. Figure 6d is GIM TEC, which is derived from CODE. It can be seen that the IRI-2016 TEC values are significantly underestimated in comparison to CODE TEC between 70 • E to 140 • E and 15 • N to 35 • N. As can be seen in Figure 6b-d, the assimilated TEC has better agreements with CODE TEC than with the IRI-2016 TEC. To a certain extent, it also indicates the improvement of the assimilation model when the GNSS-observed TEC data is assimilated into the IRI-2016 model on the disturbed day.   Figure 7a presents absolute differences between the assimilated TEC and the CODE TEC. Figure 7b presents absolute differences between the IRI-2016 TEC and the CODE TEC. Figure 7c presents assimilated TEC distribution versus the CODE TEC. Figure  7d presents IRI-2016 TEC distribution versus the CODE TEC. Figure 7a,b also present the residual histogram distributions for the data assimilation TEC and IRI-2016 TEC, respectively. It also can be seen that the residual distribution before data assimilation is obviously deviated from the center and has a large deviation from CODE TEC on the disturbed day. The residuals mean and the RMSE between the IRI-2016 TEC and the CODE TEC are −5.56 TECU and 7.53 TECU, respectively. The residual distribution after data assimilation is also closer to an unbiased Gaussian distribution. The residuals mean and the RMSE between the assimilated TEC and the CODE TEC decrease to −2.27 TECU and 3.37 TECU, respectively. In Figure 7c,d, correlation coefficients of 0.98 and 0.90 are obtained for assimilated TEC and IRI-2016 TEC, respectively. The assimilation results are also in good agreement with the CODE TEC on the disturbed day. The above research indicates that data assimilation can effectively fuse the observed data based on GNSS into the background model and obtain more reasonable and accurate results of ionospheric TEC.   Figure 7a presents absolute differences between the assimilated TEC and the CODE TEC. Figure 7b presents absolute differences between the IRI-2016 TEC and the CODE TEC. Figure 7c presents assimilated TEC distribution versus the CODE TEC. Figure 7d presents IRI-2016 TEC distribution versus the CODE TEC. Figure 7a,b also present the residual histogram distributions for the data assimilation TEC and IRI-2016 TEC, respectively. It also can be seen that the residual distribution before data assimilation is obviously deviated from the center and has a large deviation from CODE TEC on the disturbed day. The residuals mean and the RMSE between the IRI-2016 TEC and the CODE TEC are −5.56 TECU and 7.53 TECU, respectively. The residual distribution after data assimilation is also closer to an unbiased Gaussian distribution. The residuals mean and the RMSE between the assimilated TEC and the CODE TEC decrease to −2.27 TECU and 3.37 TECU, respectively. In Figure 7c,d, correlation coefficients of 0.98 and 0.90 are obtained for assimilated TEC and IRI-2016 TEC, respectively. The assimilation results are also in good agreement with the CODE TEC on the disturbed day. The above research indicates that data assimilation can effectively fuse the observed data based on GNSS into the background model and obtain more reasonable and accurate results of ionospheric TEC. Figure 8 shows the RMSE of the 3D of three modes with each epoch at different stations on the disturbed day. It can be seen that the PPP/Assimilated and PPP/CODE solutions also present a good performance as well as PPP/IF for all stations. The results also indicate that ionospheric constraint of PPP/Assimilated can effectively perform single-frequency PPP positioning under disturbed ionospheric conditions. Similar to Figure 5d, the performances of all PPP methods are generally worse in Figure 8d. In addition, this larger error occurs about 01 UT, and 13 UT to 16 UT, when ionospheric irregularity affects GNSS observations more effectively. Jiang et al. (2020) have demonstrated that ionospheric irregularities have occurred in low-latitude Southeast Asia during about 13 UT on 8 September 2017. In Figure 2, these two geomagnetic storms mainly occur at about 01 UT and 14 UT on 8 September 2017. The shock wave of the coronal mass ejection (CME) and the CME's material arrival at the earth is regarded as the main causes to induce these two geomagnetic storms [37]. The effect of ionospheric irregularities on GNSS phase observations can lead to amplitude and phase scintillation, which degrades the PPP solution. Therefore, scintillations due to ionospheric irregularities may be the main reason for the poor results of PPP/IF, PPP/CODE, and PPP/Assimilated at this time. It is interesting to find that even removing the first-order ionospheric delay in PPP/IF is not sufficient to mitigate the effects of scintillation on GNSS signals.  Figure 8 shows the RMSE of the 3D of three modes with each epoch at different stations on the disturbed day. It can be seen that the PPP/Assimilated and PPP/CODE solutions also present a good performance as well as PPP/IF for all stations. The results also indicate that ionospheric constraint of PPP/Assimilated can effectively perform single-frequency PPP positioning under disturbed ionospheric conditions. Similar to Figure 5d, the performances of all PPP methods are generally worse in Figure 8d. In addition, this larger error occurs about 01 UT, and 13 UT to 16 UT, when ionospheric irregularity affects GNSS observations more effectively. Jiang et al. (2020) have demonstrated that ionospheric irregularities have occurred in low-latitude Southeast Asia during about 13 UT on 8 September 2017. In Figure 2, these two geomagnetic storms mainly occur at about 01 UT and 14 UT on 8 September 2017. The shock wave of the coronal mass ejection (CME) and the CME's material arrival at the earth is regarded as the main causes to induce these two geomagnetic storms [37]. The effect of ionospheric irregularities on GNSS phase observations can lead to amplitude and phase scintillation, which degrades the PPP solution. Therefore, scintillations due to ionospheric irregularities may be the main reason for the poor results of PPP/IF, PPP/CODE, and PPP/Assimilated at this time. It is interesting to find that even removing the first-order ionospheric delay in PPP/IF is not sufficient to mitigate the effects of scintillation on GNSS signals.

Discussion
To verify the reliability of the assimilated method, Figure 9 shows comparisons of the TEC obtained from the GNSS observation (black) from the IRI-2016 model (blue), and from the assimilated method of the improved ensemble Kalman filter (red) at different stations during a 7-day period from 5 September 2017 to 11 September 2017. It is worth

Discussion
To verify the reliability of the assimilated method, Figure 9 shows comparisons of the TEC obtained from the GNSS observation (black) from the IRI-2016 model (blue), and from the assimilated method of the improved ensemble Kalman filter (red) at different stations during a 7-day period from 5 September 2017 to 11 September 2017. It is worth noting that these four GNSS stations are not involved in ionospheric assimilation modeling. The purpose is to use these four GNSS stations for verification of external coincidence accuracy, as the IPP near these four stations during the assimilation process will improve the modeling quality of this area. The BDS GEO satellite has a fixed IPP and can continuously observe the variation of long-term ionospheric TEC at a certain point. In this paper, we use BDS GEO satellite to analyze the accuracy of the ionospheric assimilation experiment. In Figure 9, the assimilated TEC has better agreements with observed TEC of BDS GEO than with the IRI-2016 TEC for all GNSS stations during the 7-day period. In general, the values of ionospheric TEC increase with decreasing latitude. It is consistent with the variation of ionospheric TEC. The excellent agreement between the assimilated model and the BDE GEO observations becomes even more impressive when the observed TEC values are compared with the IRI-2016 TEC values. The specific information is shown in Figure 10. Figure 10 shows assimilated and IRI-2016 TEC distributions versus the GNSS-observed TEC of different GNSS stations for a 7day period from 5 September 2017 to 11 September 2017. In Figure 10a, the correlation coefficient and the RMSE between the IRI-2016 TEC and the observed TEC of the GAMG station is 0.68 and 8.98 TECU, respectively. After data assimilation, the correlation coefficient increases to 0.92, and the RMSE decreases to 5.50 TECU. In Figure 10b, the correlation coefficient and the RMSE between the IRI-2016 TEC and the observed TEC of the JFNG station is 0.82 and 8.33 TECU, respectively. After data assimilation, the correlation coefficient increases to 0.95, and the RMSE decreases to 4.59 TECU. In Figure 10c, the correlation The specific information is shown in Figure 10. Figure 10 shows assimilated and IRI-2016 TEC distributions versus the GNSS-observed TEC of different GNSS stations for a 7-day period from 5 September 2017 to 11 September 2017. In Figure 10a, the correlation coefficient and the RMSE between the IRI-2016 TEC and the observed TEC of the GAMG station is 0.68 and 8.98 TECU, respectively. After data assimilation, the correlation coefficient increases to 0.92, and the RMSE decreases to 5.50 TECU. In Figure 10b, the correlation coefficient and the RMSE between the IRI-2016 TEC and the observed TEC of the JFNG station is 0.82 and 8.33 TECU, respectively. After data assimilation, the correlation coefficient increases to 0.95, and the RMSE decreases to 4.59 TECU. In Figure 10c, the correlation coefficient and the RMSE between the IRI-2016 TEC and the observed TEC of the LHAZ station is 0.84 and 10.64 TECU, respectively. After data assimilation, the correlation coefficient increases to 0.97, and the RMSE decreases to 5.31 TECU. In Figure 10d, the correlation coefficient and the RMSE between the IRI-2016 TEC and the observed TEC of the CMUM station is 0.93 and 12.65 TECU, respectively. After data assimilation, the correlation coefficient increases to 0.97, and the RMSE decreases to 8.21 TECU. The correlation coefficient after assimilation is increased by about 17% compared with that before assimilation, and the RMSE after assimilation is decreased by about 42% compared with that before assimilation. They demonstrate a generally good agreement between the BDS GEO TEC and the assimilated TEC. The ionospheric data assimilation can be effectively performed by the model. For further evaluation of the accuracy of the assimilated results, Figure 11 shows t RMSE of 3D with each day at different stations for days 248 to 254, 2017. In general, Figu  11 presents the gradual increase in the RMSE of the 3D coordinate for all selected statio as the latitude decreases. For stations GAMG and CMUM, the positioning accuracy PPP/Assimilated method is generally lower than that of PPP/CODE method. As show by the distribution of GNSS stations in Figure 1, we believe that this result may be due the fact that there are few stations around station GAMG and station CMUM participati in ionospheric data assimilation, resulting in low accuracy in this area. Thus, the positio ing accuracy PPP/Assimilated method is reduced. However, the positioning accuracy PPP/Assimilated method is generally higher than that of PPP/CODE method for statio JFNG and LHAZ. It is interesting to find that there are more GNSS stations around stati JFNG and station LHAZ participating in ionospheric assimilation modeling. It presen that the ionospheric model after assimilation in this area is more accurate. Figure 11 ind cates that the accuracy of 3D coordinates is improved when using the assimilated TEC the areas where there are more stations participating in assimilation. Otherwise, assim For further evaluation of the accuracy of the assimilated results, Figure 11 shows the RMSE of 3D with each day at different stations for days 248 to 254, 2017. In general, Figure 11 presents the gradual increase in the RMSE of the 3D coordinate for all selected stations as the latitude decreases. For stations GAMG and CMUM, the positioning accuracy of PPP/Assimilated method is generally lower than that of PPP/CODE method. As shown by the distribution of GNSS stations in Figure 1, we believe that this result may be due to the fact that there are few stations around station GAMG and station CMUM participating in ionospheric data assimilation, resulting in low accuracy in this area. Thus, the positioning accuracy PPP/Assimilated method is reduced. However, the positioning accuracy of PPP/Assimilated method is generally higher than that of PPP/CODE method for stations JFNG and LHAZ. It is interesting to find that there are more GNSS stations around station JFNG and station LHAZ participating in ionospheric assimilation modeling. It presents that the ionospheric model after assimilation in this area is more accurate. Figure 11 indicates that the accuracy of 3D coordinates is improved when using the assimilated TEC in the areas where there are more stations participating in assimilation. Otherwise, assimilated TEC demonstrates a reduced accuracy. The good results of the PPP/IF method indicate that almost all of these errors are efficiently mitigated. However, the RMSE of the PPP/IF method is 0.93 m, and the performance of three PPP methods is worse on disturbed day 251, 2017 in Figure 11d. As analyzed in Figure 8 above, scintillations due to ionospheric irregularities may be the main reason for the poor results of PPP/IF, PPP/CODE, and PPP/Assimilated on the disturbed day. As mentioned in Jiang et al. [38], significant ionospheric scintillation occurs in low-latitude Southeast Asia on the disturbed day. Generally, the experimental results indicate that the PPP/Assimilated method can improve the GNSS positioning accuracy more effectively in comparison to the PPP/CODE.

Conclusions
A new assimilated algorithm is successfully implemented and evaluated over China and surrounding areas. The performance of the LEnKF method is evaluated using the CODE and BDS GEO for the estimation TEC. Additionally, an application that assesses the performance of the assimilated algorithm is used to correct the kinematic single-frequency PPP. The results of statistical analysis show that correlation coefficient after assimilation is increased by about 17% compared with that before assimilation.
The RMSE after assimilation is decreased by about 42% compared with that before assimilation. Furthermore, the ionospheric constraint of PPP/Assimilated can effectively improve the accuracy of single-frequency PPP. On the disturbed day, we find that even removing the first-order ionospheric delay in PPP/IF is not sufficient to mitigate the effects of scintillation on GNSS signals. Scintillations due to ionospheric irregularities may be the main reason for the poor results of PPP/IF, PPP/CODE, and PPP/Assimilated methods. Therefore, the LEnKF method can be considered as a useful tool of data assimilation for

Conclusions
A new assimilated algorithm is successfully implemented and evaluated over China and surrounding areas. The performance of the LEnKF method is evaluated using the CODE and BDS GEO for the estimation TEC. Additionally, an application that assesses the performance of the assimilated algorithm is used to correct the kinematic single-frequency PPP. The results of statistical analysis show that correlation coefficient after assimilation is increased by about 17% compared with that before assimilation.
The RMSE after assimilation is decreased by about 42% compared with that before assimilation. Furthermore, the ionospheric constraint of PPP/Assimilated can effectively improve the accuracy of single-frequency PPP. On the disturbed day, we find that even removing the first-order ionospheric delay in PPP/IF is not sufficient to mitigate the effects of scintillation on GNSS signals. Scintillations due to ionospheric irregularities may be the main reason for the poor results of PPP/IF, PPP/CODE, and PPP/Assimilated methods. Therefore, the LEnKF method can be considered as a useful tool of data assimilation for applications of space environment monitoring and navigation positioning. In the future, we will further study the multi-scale regional ionospheric assimilation model. Meanwhile, it will be further evaluated whether the multi-scale ionospheric assimilation model can improve the accuracy of single-frequency PPP, especially under conditions of ionospheric anomalies.