An improved method using adaptive smoothing for GNSS tomographic imaging of ionosphere

Global navigation satellite system (GNSS) is a well-established sensors in the recent ionosphere research. By comparing with classical meteorological equipments, the GNSS application can obtain more reliable and precious ionospheric total electron content (TEC) result. However, the most used GNSS ionospheric tomography technique is sensitive to a priori information due to the sparse and non-uniform distribution of GNSS stations. In this paper, we propose an improved method based on adaptive Laplacian smoothing and algebraic reconstruction technique (ALS-ART). Compared with traditional constant constraints, this method is less dependent on a priori information and adaptive smoothing constraints is closer to the actual situation. Tomography experiments using simulated data show that reconstruction accuracy of ionospheric electron density using ALS-ART method is significantly improved. We also use the method to do the analysis of real observation data and compare the tomography results with ionosonde observation data. The results demonstrate the superiority and reliability of the proposed method compared to traditional constant constraints method which will further improve the capability of obtaining precious ionosphere TEC by using GNSS.


Introduction
The ionosphere is an ionized region of the Earth's atmosphere, which is generally accepted that begins at 50 km and ends at 1000 km approximately from the Earth surface. In order to convert accurate integral measurements into 2-D or 3-D structures, a technique called tomography has been invented. In the field of global navigation satellite system (GNSS) atmosphere, the principle of computerized ionospheric tomography (CIT) becomes applicable with the increasing number of GNSS satellites and the build-up of ground-based GNSS stations in the 1990s [1][2][3][4]. Since then, a variety of CIT approaches based on GNSS observation data have been developed for the accurate reconstruction of ionospheric electron density (IED) distribution in the upper atmosphere [5][6][7][8][9][10][11][12][13][14]. An overview about the direction and challenges of an area at the forefront of CIT research is provided by Bust and Mitchell [15]. CIT  inverse problem which is ill-posed with the limited number of viewing angles of the raypath measurement. To get the uniqueness of the CIT solution, there are many approaches which are used to solve the reconstruction sparseness.
Smoothing constraints approach can effectively constrain near voxels with no raypath information. Garcia and Crespon [16] applied local basis parametrisation of electron density which constrained by NeQuick ionosphere model and its spatial gradients. Hobiger et al. [17] used constrained simultaneous algebraic reconstruction technique for tomography reconstruction of electron density distribution and found that the convergence speed would be significantly higher than that of classical method. Lee and Kamalabadi [18] regularized the inverse problem by incorporating neighborhood smoothness and continuity constraints applicable to general ionospheric conditions. Wen et al. [19,20] proposed constrained algebraic reconstruction technique for tomographic reconstruction and demonstrated the reliability and superiority of this method. Nesterov and Kunitsyn [21] developed minimal Sobolev's norm constraint to smooth the solution. Yao et al. [22] incorporated Tikhonov and total variation regularization constraint to overcome near voxels with no raypath information. Panicciari et al. [23] permited sparsity in the inversion coefficient by using minimization constraint and wavelet basis functions. Norberg et al. [24,25] used Bayesian statistical inversion with prior distribution given by its mean and covariance with ionosonde measurements and Gaussian Markov random fields as a sparse matrix for the numerical computations, and employed the required additional information with Bayesian which could be given with prior probability distributions, respectively. Wang et al. [26,27] used model function in regularization constraint to balance the weights between observations and electron density field, and model function in the modified L-curve method to balance a prior information and real measurements, respectively. To a certain extent, these existing studies have overcome the ill-posedness of ionospheric electron density reconstruction based on constraints and priori information, especially in vertical structure of the ionosphere. These studies have shown that the results of CIT could strongly depend on constraints and prior information.
However, if the reconstruction relies excessively on a priori information in the CIT process, the tomographic results are similar to the prior information, and the CIT technology will lose its inherent value. To reduce the dependence on the prior information and reconstruct a higher precision ionospheric vertical structure, we propose a novel approach using adaptive Laplacian smoothing and algebraic reconstruction technique (ALS-ART) for three-dimensional GNSS ionospheric tomography. In the following, a more detailed discussion of CIT constraints is provided. Therefore, Section 2 describes the models for GNSS measurement and the principles of ionospheric tomography. Section 3 introduces how the adaptive smoothing is used for ionospheric electron density (IED). In Section 4, the proposed approach is applied to simulated and real tomography and this section also discusses the problems of CIT. Section 5 concludes the major findings.

Ionospheric total electron content extraction
Ionospheric total electron content (TEC) is the integrated ionospheric electron density along the ray path between GNSS satellites and receivers. The TEC is commonly extracted from dual frequency GNSS receivers through methods of carrier phase smoothed pseudorange [28,29] and uncombined precise point positioning (UPPP) [30,31]. The TEC estimated by the method of carrier phase smoothed pseudorange is affected by code-delay multi-path and length of every continuous arc [32,33]. Zhang et al. [30] proposed the UPPP method for estimating ionospheric observables, and they found that it effectively reduces leveling errors. The detailed introduction of the method can be found in Zhang et al. [30]. In this paper, we use the method to extract ionospheric TEC. The observation equations of the GNSS pseudorange and carrier phase can be generally expressed as where P s r;f and L s r;f are pseudorange and carrier phase observation from receiver r to satellite s on frequency f(f = 1,2) in length, respectively; ρ is the geometry distance between the receiver and satellite; c is the speed of light; δt γ and δt s are the receiver and satellite clock offset, respectively; T trop and I ion are the tropospheric and ionospheric delay; λ is the wavelength; b is the frequency dependent code bias and d is the frequency dependent phase bias; N is the integer phase ambiguity in cycles; ε is the observation noise.

Tomography model
The slant TEC data are used as CIT projections in a region [34], which along the ray path of the GNSS signal between a satellite and a ground receiver is defined as the integrated value of the IED and is modeled as where STEC(t) is the slant TEC; Ne(r,t) is the electron density at the time t; I and J are the respective total number of receivers and satellites; and i and j are the positions of the ith ground receiver and jth satellite, respectively. In a CIT approach, the ionosphere is divided into a grid in which each voxel has the same size or is not equidistant from other voxels. The linear integral in Eq (2) includes the slant TEC contribution along the entire ray path from the satellite to the receiver [18][19][20][21][22]. A parametric representation of Eq (2) can be written as where n and a denote a sampling point and the weights for numerical integration at the sampling points, respectively; N is the total number of voxels in the ionosphere grid; and ε is the measurement and model error.
The formulation in Eq (3) can be expressed as a discrete mathematical problem written as where Y 2 R M×1 is the slant TEC from GNSS observation values, M is the total number of slant TEC measurements, A 2 R M×N is the observation matrix that corresponds to the grid, X 2 R N×1 is the IED at each voxel, and ε denotes the noise.

Classic method and problem
Algebraic reconstruction technique (ART) is a row-action technique which is a classical iterative algorithm. It sets a priori information for each voxel before the iteration begins in the tomographic region. The kth iteration of ART algorithm computes the difference between Y and Y k which is obtain by using the current estimated value X k in Eq (4) [35][36][37]. A correction derived from the difference is distributed over X k to get X k+1 . After numbers of iterations, the tomography result converges to a solution of Eq (4). For the kth iteration, the ART algorithm can be written as where A i is the ith row of matrix A. k is the number of iterations. γ k is relaxation parameters which is usually chosen to be the same and confined to the interval (0,2). Generally, the number of rays in tomography is larger than the number of voxels, so the ray geometry defines a non-singular matrix, which permits reconstruction of the detected media [17]. However, this is not usually the case in geoscience applications. For the ionospheric tomography, there are near voxels with no raypath information and the final values are often the same as a priori information. For this problem, additional constraints or a priori information is generally used for ionospheric tomography [17,19]. The usual way to constrain the inverse problem is through damping and smoothing once it is linearized in Eq (4). A voxel without any ray passing through will extract information from their near voxels, while an over-determined voxel will determine their values through the data and contribute to the determination of their near voxels. These constraints can be written as where H, V, B are the horizontal constraints, vertical constraints and boundary conditions. According to Eqs (4) and (6), we can get The Eq (7) can be expressed as where S includes A, H, V and B. From Eqs (5) and (8), the algorithm combined by smoothing and ART can be given as The application of smoothing constraints indirectly increases the amount of observation information, which is equivalent to adding virtual observations to the voxels without any observation information, overcoming the problem of excessive dependence on a prior information.

Constant Laplacian smoothing
Before introducing the new tomography algorithm, we will present Laplacian smoothing which is the basis of the new algorithm. Laplacian smoothing is an algorithm that can smooth polygon meshes based on temporal or spatial continuity. The smoothing operation is as follows: The value x 0 of the voxel is replaced by the average of its neighboring voxels. For tomography rectangular grid, we only consider the nearest voxels. Because reconstruction region has strict boundaries, there are three cases when we use Laplacian smoothing, as shown in Fig 1. In Fig 1, x 1 , x 2 , x 3 and x 4 are the nearest voxels to x 0 in a CIT region. We use Laplacian smoothing method to constrain them. It can be expressed as in Fig 1(a), in Fig 1(b) and in Fig 1(c).
where q is smoothing factor whose initial values are 4, 3 and 2 for Fig 1(a), 1(b) and 1(c), respectively. The three equations are used to construct each voxel as constant smoothing constraints in the horizontal and vertical directions, respectively. However, constant smoothing factors may be inaccurate. Thus, we use an adaptive Laplacian smoothing method to reconstruct IED distribution.

Adaptive Laplacian smoothing
If unreasonable constraints are imposed, the accuracy of the IED reconstruction will be significantly reduced. To make the constraints closer to the actual situation, we develop an adaptive smoothing algorithm by interacting between constraints and reconstruction results. Firstly, we use the Laplacian smoothing to impose the initial horizontal and vertical constraints and use the ART algorithm to inverse the initial solutions. Then, we use the initial solutions to estimate the smoothing factors q with Eqs (10) to (12) and establish the new horizontal and vertical constraints. The new ionospheric tomography approach ALS-ART is performed as follows.

Calculate IED values by Eq
where m is 4, 3 and 2 from Eqs (10) to (12), respectively. x h is a threshold to prevent updating the smoothing factor from inaccurate results and is a half of the maximum IED in the results. We set a scale factor to update the x h until it is less than triple root mean square error of IED. Repeat steps (2)-(3) until the difference between the IED values are smaller than the empirical 0.3 TECu.

Experimental results
To

Simulation
In this section, the real GNSS geometry used to construct matrix A in Eq (4) is used to calculate the slant TEC values for the simulated ionosphere. The simulated IED distributions are from  Figs 3 and 4, we can see that the accuracy of reconstructing IED using CLS-ART is obviously lower than that of reconstructing IED using ALS-ART between 100 km and 400 km. Table 1 provides the reconstruction results of the error statistics for the two methods based on various longitudes. The maximum absolute error, average absolute error and RMS error of IED with the ALS-ART algorithm are smaller than those of IED with the CLS-ART algorithm at various longitudes.   reconstructing IED using CLS-ART is obviously lower than that of reconstructing IED using ALS-ART at 150 km, 250 km, 350 km and 450 km of altitude. Table 2 represents the error analysis of IED reconstruction by using the two methods at various heights, and the results also presents the reconstruction results of the error statistics for the two methods based on various altitudes. The maximum absolute error, average absolute error and RMS error of IED with the ALS-ART algorithm are also smaller than those of IED with the CLS-ART algorithm at various longitudes. These results indicate that the CIT accuracy of the ALS-ART algorithm is generally superior to that of the CLS-ART algorithm. In addition, the maximum absolute error, average absolute error and RMS error of IED with the ALS-ART algorithm based on total voxels are approximately 3.51×10 10 el/m 3 , 0.25×10 10 el/m 3 and 0.36×10 10 el/m 3 , respectively. The maximum absolute error, average absolute error and RMS error of IED with the CLS-ART algorithm based on total voxels are approximately 6.73×10 10 el/m 3 , 0.43×10 10 el/m 3 and 0.52×10 10 el/m 3 , respectively. We believe that the ALS-ART algorithm has the validity and superiority. Additionally, we think that the ALS-ART algorithm can derive more reliable values and obtain results closer to the true values compared with the CLS-ART algorithm. We use observed slant TEC data from the IGS in Europe with a vertical profile derived from the two ionosonde stations. We use the real measurements from 22 August 2018, when the ionosphere was under calm conditions. Fig 7 demonstrates the hourly variation of the IED at a fixed longitude meridian of 10˚E on 22 August 2018. This figure shows that the IED characteristics of the different latitudes are quite different, and the values of IED in the north is smaller than those of IED in the south as a whole. The characteristics indicate a strong correlation between IED and latitude. Comparing all the sub-graphs in Fig 7, it can be seen that the peak height of the IED gradually decreases in the period of 2:00 UT to 10:00 UT. Then, it gradually increases in the next time period. The peak height of the IED rises to 350 km at 22:00 UT, which presents the characteristics of the vertical variation of the IED this day.

Real observation data
To verify the reliability of the CIT results, we use the ALS-ART algorithm to inverse IEDs which are compared with ionosonde data.  The maximum absolute error of IED with the CLS-ART algorithm based on Dourbes and Pruhonice are approximately 1.  This result shows that the inversion accuracy of the ALS-ART algorithm is better than that of the CLS-ART algorithm. Fig 9 shows a comparison of the F2 layer peak electron density (NmF2) obtained from the two algorithms, as well as the observation data from the ionosonde stations at hourly time intervals on 22 August 2018. The peak in the electron density of the F2 layer obtained by two algorithms is generally in good agreement with that based on ionosonde measurements, but the peak obtained by the ALS-ART algorithm is generally closer to the ionosonde measurements than that based on CLS-ART algorithm.
For further validation, we give several examples of CIT results for regional IED distribution during the geomagnetic storm on 26 August 2018. The results are represent in Fig 9. Fig 10 demonstrates the hourly variation of the IED at a fixed longitude meridian of 10˚E on this day. This figure also shows that the IED characteristics of the different latitudes are quite different, and the values of IED in the north is smaller than those of IED in the south as a whole. The peak height of the IED gradually decreases to 250 km in the period of 2:00 UT to 10:00 UT. Then, it gradually increases in the next time period. The peak height of the IED rises to 350 km at 22:00 UT, which also presents the characteristics of the vertical variation of the IED this day.   of the IED under disturbed ionospheric conditions. An improvement in the NmF2 values is very obvious, but the peak height hmF2 values are apparently not well reconstructed. In Fig 11 (a) and 11(b), the peak height based on the ALS-ART algorithm is slightly higher than that of the ionosonde data. Fig 12 shows the characteristics of the reconstructed F2 layer are validated against those of the ionosonde measurements. It is clear that the ALS-ART algorithm can improve the estimated NmF2 values with respect to the CLS-ART algorithm. Due to GNSS observation noise and station geometry limitations, the vertical resolution of the CIT results still needs to be improved. To evaluate reliability of the proposed method, we show the F2 layer peak density (NmF2) using the two methods with one week data. As shown in Fig 13, the characteristics of the reconstructed NmF2 are validated against those of the ionosonde measurements. Generally, the AS-ART method can improve the estimated NmF2 values with respect to the CS-ART method.
Reconstruction of tomographic images based on GNSS is a mathematical inversion problem. Iterative algorithms are commonly applied, but it is not well understood how good the solutions are. Depending on a prior information, the iteration may lead to different results which all satisfy the original measurements equally well [19]. The quality of the reconstruction results depends on how suitable the prior information are. The solution of smoothing constraints is well understood. It presents the most probable values of the unknowns by conditional constraints once the voxels without any observation information. One drawback of the smoothing operator is that, for applications in which IED spatial distribution contain discontinuities, the CLS-ART algorithm could present incorrect results. In the case of steep gradients, the adaptive smoothing operator should be used in order to allow for more variation between voxels. In this study, we have proposed a novel approach ALS-ART and tested the results from this algorithm inversion of CIT against independent ionosonde measurements. In general, the inversion accuracy of the ALS-ART algorithm is better than that of the CLS-ART algorithm. An improvement in the NmF2 values is very obvious. The inversion results of the method are different from the ionosonde measurements in the vertical direction because the grid should be unrefined in the vertical direction.

Conclusion
In this study, we implement the novel approach based on adaptive smoothing and algebraic reconstruction technique (ALS-ART) and conduct a series of ionospheric tomography experiments to investigate the accuracy of inversion results using the algorithm. The algorithm imposes adaptive constraints according to the smoothness of near voxels. To some extent, it overcomes the disadvantages of the CLS-ART algorithm. Experiments involving simulated values and real measurements are performed to validate the feasibility and superiority of the proposed approach compared to CLS-ART algorithm. The ALS-ART algorithm have been used to reconstruct IED under geomagnetic calm and storm conditions. The IED profiles inversed from ALS-ART match better with ionosondes than those from CLS-ART. Although the method improves the ionospheric tomography results, it cannot completely improve the vertical resolution of tomography. However, it is necessary to further add the quantity of observation information.