Three-dimensional ionospheric tomography based on compressed sensing

The insufficient number of low-elevation observations is a limitation of the three-dimensional ionospheric computer tomography (CT) based on the global navigation satellite system (GNSS). To solve this problem, accurate prior information on the regional ionosphere must be obtained. However, it is difficult to explicitly and accurately express prior ionospheric information. This study uses compressed sensing (CS) for ionospheric tomography for the first time. Specifically, the electron density obtained from the international reference ionosphere is used to build a dictionary to fully integrate the prior information into the dictionary. Then, the electron density is reconstructed by using the compressive sampling matching pursuit method. Subsequently, the GNSS data of China (Region I) and Europe (Region II) were utilized to validate this proposed method, and the results are compared with ionosonde observations. The mean and standard deviation (SD) of the difference with respect to the ionosonde result are 41 and 22 km, respectively. The mean and SD of relative deviation were 16% and 9%, respectively. In Region II, the mean and SD of the deviation between the reversed peak electron density and the result of the ionosonde were 1.9 × 1010 m−3 and 8.1 × 1010 m−3, respectively. The mean and SD of the relative deviation were 3% and 13%, respectively. The mean and SD of the peak height deviation were 33 and 19 km, and the mean and SD of the relative deviation were 11% and 7%. The electron density distribution and variation in these two regions showed a local time dependence, and the horizontal gradient of the electron density in the latitude was greater than that in the longitude. Moreover, CT by CS is efficient, taking about 6 s per inversion based on an desktop computer with 16 GB RAM and Intel (R) Core (TM) i7-8700 CPU.


Introduction
Using the total electron content (TEC) in ionospheric tomography was first suggested by Austen et al. (1986). Andreeva et al. (1990) successfully implemented experimental radiotomographic reconstructions. Xu et al. (1995) carried out two-dimensional ionospheric tomography experiment at low latitudes in East Asia to detect ionospheric anomalies. Since the late 1990s, the global navigation satellite system (GNSS) has been widely used in ionospheric research due to the high temporal and spatial resolutions as well as the low cost, which led to the era of the three-dimensional ionospheric computer tomography (CT) (Rius et al. 1997). Thus, the GNSS measurements became the main data source in three-dimensional ionospheric tomography. However, threedimensional ionospheric tomography has certain limitations, especially for the insufficiency of horizontal observations (Yeh and Raymund 1991). Numerical inversion algorithms have been developed to overcome this problem. Generally, ionospheric tomography algorithms can be classified into two categories. The first is pixel-based algorithms, including additive algebraic reconstruction (Austen et al. 1988;Xu et al. 2005), multiplicative algebraic reconstruction (Raymund et al. 1990), simultaneous iterative reconstruction (Pryse et al. 1993), and the method of constrained least squares (Saito et al. 2017). To improve the iterative algorithm, Yao et al. (2015) imposed a priori constraints by increasing the virtual observations, Yao et al. (2018) applied the side rays to the inversion, and Zheng et al. (2018) proposed an automatic search technology of relaxation factor. The second is function-based algorithms, including the orthogonal function method and singular value decomposition (Fremouw et al. 1992;Raymund et al. 1994). Farzaneh and Forootan (2018) improved the empirical orthogonal function method by combining the spherical Slepian function with the orthogonal empirical function. In recent years, some new methods have been applied to CT. Harding and Milla (2013) applied compressed sensing (CS) to one-dimensional imaging of coherent backscatter from ionospheric plasma density irregularities at the magnetic equator. Panicciari et al. (2015) proposed a CT algorithm based on l 1 minimization for CT by using GNSS phase measurement. Hysell et al. (2019) evaluated compressed sensing (CS) methods in the application of two-dimensional aperturesynthesis imaging of radar backscatter from field-aligned plasma density irregularities in the ionosphere, the evaluated CS methods included basis pursuit denoising, implemented with the fast iterative shrinkage thresholding algorithm, and orthogonal matching pursuit with a wavelet basis. Sui et al. (2021) applied compressed sensing (CS) based on l 1 norm for the sparse reconstruction of 3-D regional ionospheric tomography using the differential STEC.
The theoretical framework of CS was constructed in 2006 (Candès and Tao 2006;Donoho 2006). CS shows that a signal can be reconstructed from a small number of measurements if the signal is sparse in a fixed basis or compressible. CS includes three parts: the sparse representation of the signal, the requirements for the observation matrix, and the reconstruction algorithm.
The sparse representation of the signal means that the signal only contains a few large and sparse elements, or it can be converted into a sparse vector using a matrix which is called a dictionary. Although most of the signals are not sparse, they generally follow a certain regularity. Thus, they can be sparsely represented by a dictionary. The two types of dictionaries are the complete orthogonal dictionary and the over-complete dictionary (Mallat and Zhang 1993). The complete orthogonal dictionary includes Fourier basis and discrete cosine transform dictionary. The over-complete dictionary can be further divided into fixed and learning dictionaries. Once a fixed dictionary is confirmed, it will not change, such as the Gabor dictionary (Bergeaud and Mallat 1998), Gaussian dictionary (Qian and Chen 1994) and cascade dictionary (Elad 2010). A learning dictionary is obtained through learning or training. Typical learning methods include the method of optimal directions algorithm (Elad and Aharon 2006), the recursive least squares algorithm (Skretting and Engan 2010), KSVD (Aharon et al. 2006), and the double sparse method (Rubinstein et al. 2010). Candès and Tao (2005) proposed that the observation matrix should satisfy the restricted isometry property (RIP). Baraniuk et al. (2008) found that the observation matrix was not related to the atoms in the dictionary, which was equivalent to RIP. Candès and Plan (2011) suggested that the observation matrix selects sensing vectors independently at random from a probability distribution F, which includes all standard models and a framework for new measurement strategies, and the probability distribution F obeys a simple incoherence property and an isotropy property. However, it is difficult to determine whether the observation matrix satisfies RIP. Candès and Plan (2011) also suggested that RIP was too strict and conservative. Zhang (2013) demonstrated that, based on prior knowledge, the signal could be accurately reconstructed even if RIP was not satisfied. Moreover, Adcock et al. (2017) studied a situation where the observation matrix was highly correlated with the dictionary, but the signal was accurately reconstructed and proposed a framework that included the concepts of asymptotic sparsity, asymptotic incoherence, multisampling, and infinite-dimensional CS. Candès and Tao (2006) showed that RIP could most likely be satisfied when the Gaussian matrix, Fourier matrix, and binary matrix were used as the observation matrices. However, these matrices were random matrices. In some cases, it is hard to observe randomly, such as in ionospheric tomography. Therefore, it is one of the main research directions of CS to construct deterministic observation matrices (Elad 2007).
CS reconstruction algorithms can be divided into three categories. The first one is convex relaxation algorithms, in which the non-convex sparsity 0 -norm was replaced with 1 , which has the convex property, thereby obtaining a convex optimization problem that is easier to solve. Convex optimization algorithms include constrained optimization (Figueiredo et al. 2007) and approximate operator-based strategies (Parikh and Boyd 2014). The second one is the greedy algorithms, including matching pursuit (Bergeaud and Mallat 1998), orthogonal matching pursuit (Tropp 2004), the iterative hard threshold method (Blumensath and Davies 2009), and the compressive sampling matching pursuit (CoSaMP) (Needell and Tropp 2009). The third category includes focal underdetermined system solver (Rao and Kreutz-Delgado 1999) and iteratively reweighted least squares (Chartrand and Yin 2008).
In the present study, three-dimensional ionospheric tomography based on CS is proposed. The electron density obtained from IRI empirical model is used to construct a dictionary, and then, the observation matrix is optimized based on weights. Lastly, unlike the previous studies (Panicciari et al. 2015;Sui et al. 2021), the CoSaMP reconstruction method is used to invert the electron density based on STEC data in China and Europe region. The inversion results are compared with ionosonde measurements.

Method
CS is a new signal sampling and reconstruction framework that can achieve the high-precision reconstruction of the original signal with very few samples. The observation matrix, sparse representation and reconstruction algorithm are the three key components of CS. When the observation matrix satisfies certain conditions, the signal can be accurately reconstructed with a very high probability. Sparseness means that only a few elements in the signal are non-zero or a few elements have large values. Although natural signals are generally non-sparse, most of them can be converted into a sparse signal through a dictionary. In addition, there are many reconstruction algorithms in CS theory, and the CoSaMP was used in this study.

Model
A pixel-based model was used in this study. The inversion area was divided into N grids, and it is assumed that the electron density in each grid was the same. Figure 1 shows the grid division in Region I (China). The longitude range of the inversion area was 104.5°-120.5° E, the latitude range was 24.5°-40.5° N, and the altitude range was 95-2100 km. The inversion area was divided into non-uniform grids, and there were five types of grids of different sizes. The five types of grids differed only in the altitude resolution, whereas the latitude and longitude resolution was always 1° × 1°. Table 1 shows the height information of the grids.  For a GPS ray passing through the inversion area, the TEC can be expressed as: where STEC g is the TEC in the inversion area, i is the grid number, l i is the intercept of the ray in the grid, Ne i is the electron density at the center of the grid, and d and s represent the discretization error and the observation error, respectively. Based on the NeQuick model, the proportion r g of STEC g in the inversion area to the total STEC can be obtained by following the method proposed by Yao et al. (2018): where r and s are the errors of r g and STEC, respectively. By substituting (2) into (1), we arrive at: The matrix form of (3) is as follows:     where each item corresponds to the vector or matrix form of each item in (3). For example, is the column vector form of r g ⋅ STEC . From (4), it can be seen that the observation equation contains many error terms, including the discretization error caused by gridding , the scale factor error , and the observation error of STEC .

Sparse representation
A prerequisite for signal reconstruction using CS is a sparse representation of the signal. Although the regional electron density column vector is usually not sparse, it can be represented sparsely after dictionary conversion. To verify the sparsity of ionospheric electron density, an experiment was carried out based on the IRI model. The IRI model was used to obtain the long-term electron density in the inversion area, and then, the singular value decomposition was used to construct a long-term dictionary: where is the conversion result of using dictionary . The electron density at the center of the grid in the inversion area is expressed as: where is a column vector composed of the electron density at the center of the grid in the inversion area. Due to the lack of real electron density data in the entire region, the sparse representation of the ionosphere in this region was evaluated based on . Figure 2 shows the elements in the first 50 rows and first 100 columns of , where is the normalized , i.e., (i, j)= abs( (i,j)) max(abs( (∶,j))) . It can be seen that a few elements in the front rows have large absolute values, indicating that has a good sparsity and the ionospheric electron density can be sparsely represented by a, which only contains k large elements, and the other elements are set to zero. Specifically, the sparsity of a is k, and the location of all non-zero elements in a (i.e., the row number) is called the support set. Equation (7) is converted to: where is the error of a. By combining (4) and (8), we obtain: where is the observation matrix, is sparse, which makes it possible to reconstruct using CS.

Observation matrix weight
The weighting of the observation matrix is complex and is restricted by two factors: the error term and the reconstruction performance of the observation matrix. In the algorithm, the inner product and the method of the least squares were used. An appropriate weight of the observation equation can help to obtain good results. However, the weight can change the reconstruction performance of the observation matrix positively or negatively. Therefore, both factors should be considered when determining the weight.
In the CS theory, the observation matrix should satisfy certain conditions in order to reconstruct the signal with high probability and accuracy. Candès and Tao (2005) proposed that the observation matrix needed to satisfy RIP. Specifically, for k = 1, 2, 3, … , K , the isometric constant k of the matrix is defined as the minimum value that satisfies: where || || 2 is the 2 norm and x is a k column sparse vector. If 0 < k < 1 , then satisfies k-order RIP.
Traditionally, the weight is generally determined based on the covariance matrix. As can be seen from the supplementary material, the error terms include the discretization error , the scale factor error , the observation error term , and the sparse representation error term ⋅ . The discretization error can be neglected because it is on the order of 0.01 TECU. See supplementary material for detailed error analysis. In one calculation, S r ∼ N(m Sr , 2 Sr ) , Sr ≪ | | m Sr | | . r s ∼ N(0, ...) , and may not be true in one calculation. Ignore the scale factor error for its small standard deviation and the uncertainty r s , only consider the sparse representation er ror term is appropriate. However, in CS, the performance of the observation matrix directly affects the reconstruction result or even leads to an incorrect result. The performance of the observation matrix is affected by the weight. Therefore, when determining the weight, we focused on the performance of the where w is the weight, i.e., the weight of the observation is the reciprocal of the total intercept of the observation. By combining (9) with the above equation, we obtain: where is the weight matrix of w , and the off-diagonal elements in the matrix were 0.

Reconstruction algorithm
The CoSaMP algorithm was used to reconstruct the electron density. The CoSaMP algorithm was a special greedy algorithm. Theoretical performance analysis was difficult for the greedy algorithm (Blumensath et al. 2012). Davenport and Wakin (2010) reported that the theoretical conditions of the greedy algorithm were stricter than RIP, while other researchers showed that the theoretical performance of some greedy algorithms was close to and even better than that of the convex relaxation method in some specific situations (Tropp 2004;Needell and Vershynin 2010;Blumensath and Davies 2008).
CoSaMP is different from classical matching pursuit algorithms. In each iteration, the CoSaMP algorithm chooses multiple atoms in the support set based on the correlation. Then, the support set was cropped using the method of least squares. The specific steps are as follows: 1. Initialization. For = ⋅ , where and correspond to and in (13), respectively. Let = , = , the sparseness is k, and the support set is empty = � . The columns of A are normalized to obtain B. 2. The residual error is updated, and the normalized sensing matrix is correlated with the residual error. = − ⋅ , = ⋅ . 3. The set R of the 2k row numbers that are most correlated with the residual error is obtained, namely, the row R corresponding to the 2k elements with the largest absolute value in p, = ∪ . 4. b is updated using the least squares, and the support set is cropped. ( ) = ( ) −1 ⋅ + ( ) . The row corresponding to the k elements with the largest absolute value was obtained from b, = . 5. Update b. ( ) remains constant, and the other rows of b are set to 0. 6. Steps 2-5 are repeated until ry no longer changes, or the 1-norm of ry is less than the threshold, which is set to 5% of the 1-norm of y. The results are calculated. Using the final support set, the sparseness representation of the electron density a is calculated using the least squares method. Then, the electron density is obtained using (8).

Dataset
GPS observations on March 10, 2015, in China region (Region I) and on July 5, 2014, in the Europe region (Region II) are used to reconstruct the ionospheric electron based on CS. The Kp and Dst index is presented in Table 2 and Fig. 3. The STEC data of the Region I were derived from the terrestrial network by code-phase leveling, and DCBs were estimated by the spherical harmonic model (Jin et al. 2012). The STEC data of Region II can be downloaded from http:// www. gage. upc. edu/ produ cts (Rovira-Garcia et al. 2016a, 2016b. Ionosonde observations at Wuhan (30.50°, 114.40°) and Dourbes (50.10°, 4.60°) are also used to compare with the reconstructed results. Ionosonde data can be downloaded from https:// lgdc. uml. edu/ common/ DIDBF astSt ation List. Figure 4 shows the three-dimensional map of Region II, and Figs. 5 and 6 show GPS stations, inversion area, and ionosonde stations used in Region I and Region II, respectively.

Results and discussion
Figures 7 and 8 compare the variation of the electron density peak height reconstructed by CS (blue circles), generated from IRI (black line) and measured by Wuhan and Dourbes ionosondes (red dots) in Regions I and II, respectively. We can see from Fig. 7 that in Region I the peak height reconstructed by CS shows a similar trend to that generated from the IRI model and measured by Wuhan ionosonde, and the peak height reconstructed by CS is higher than that obtained from the other two ways, this may be caused by the insufficient applicability of the dictionary. As shown in Fig. 8, the peak height reconstructed by CS was higher than that calculated by IRI model and observed by Dourbes ionosonde in Region II. Besides, the trend of the peak height reconstructed by CS is consistent with that obtained from the IRI model and Dourbes ionosonde. Then, we take the Wuhan and Dourbes ionosonde measurements as the reference and calculate the error and SD of the peak height of the IRI model and CS, as shown in Table 3. In Region I and Region II, the mean error and SD of CS are larger than those of IRI model, which means the IRI model can predict the peak height more accurately.
Figures 9 and 10 compare the peak electron density reconstructed by CS (blue circles), predicted by IRI model (black line) and observed by ionosonde (red dots) in Region I and Region II. As shown in Fig. 7, in Region I, the peak electron density reconstructed by CS was very close to that observed by the ionosonde, while IRI model predicted result is much larger than that of CS and Wuhan ionosonde observation. In Region II, the CS result shows a similar trend to the ionosonde observation. Table 4 shows the mean error and SD for the peak electron density in Region I and Region II, and the ionosonde observations at Wuhan and Dourbes were taken as the reference. Obviously, in Region I, the mean error and SD of the peak electron density reconstructed by CS are smaller than those predicted by the IRI model, which indicates that CS can reconstruct the peak electron density more accurately. However, in Region II, the  Electron density sections along the longitude direction in Region II at 00:00 UTC time difference in the mean error and SD between CS and IRI models is not obvious. Figure 11 shows the comparison of the electron density profiles reconstructed by CS (blue lines), predicted by IRI model (black lines) and inversed by the SAO Explorer software based on the ionosonde observations at Wuhan (red lines) during the period from 00:00 to 23:00 UTC on March 10, 2015, in the Region I. It needs to be noted that the profile above the peak height of the ionosonde was extrapolated. As shown in Fig. 11, the electron density profiles reconstructed by CS are in good agreement with that observed by the Wuhan ionosonde, although the peak height of the electron density reconstructed by CS is generally higher than that observed by an ionosonde. IRI model predicted electron density profiles are twice the ionosonde observations. Figure 12 compares the electron density profiles reconstructed by CS (blue lines), predicted by IRI model (black lines) and inversed by the SAO Explorer software based on the ionosonde observations at Dourbes (red lines) during the period from 00:00 to 23:00 UTC on July 5, 2014, in the Region II. These three profiles are very similar. In addition, compared with Fig. 11, the reconstruction result of CS in Region II was superior to that in Region I, especially before 11:00 UTC. Figure 13 shows the electron density slices at different latitudes at 00:00 UTC on March 10, 2015, in Region I by CS. All slices show enhanced electron density in the region between 200 and 400 km, and the electron density of this area decreases with the increase of the latitude. In addition, the electron density also enhances with the increase of the longitude. This is because Region I is in the morning, and the solar radiation was stronger in large-longitude and lowlatitude areas. Figure 14 shows the electron density slices at different longitudes at 00:00 UTC in Region I by CS. Similar to Fig. 13, enhanced electron density was observed between 200 and 400 km. Also, the electron density in this region decreased with the increase of the latitude and enhanced with the increase of the longitude. Figure 15 shows the electron density slices at different altitudes in Region I at 00:00 UTC by CS. We can see from Fig. 14 that the maximum electron density in the slice at 300 km appears at 25° latitude and 120° longitude, which was consistent with the findings of Figs. 12 and 13. Figures 16, 17 and 18 present the time variation of the electron density slices at 32° latitude, 112° longitude, and 300 km altitude in Region I by CS, respectively. All these three figures show that the electron density increases with time and peaks at 08 UTC; after that, it decreases rapidly. It can be seen from Fig. 16 that the enhanced electron density area from 200 to 400 km moved from a high-longitude area to a low-longitude area with time and then disappeared at 16:00 UTC. As presented in Fig. 17, the enhanced electron density area expanded from low-latitude areas to highlatitude areas and then vanished at 16:00. The same phenomenon is shown in Fig. 18 as well. It can be seen from Figs. 16 and 17 that the changes in electron density along the longitude direction were smaller than that along the latitude direction. In addition, the enhanced electron density area between 200 and 400 km, which appeared at 0:00 UTC, and the center of the spheroid moved to 112°longitude and below 25° latitude at around 8:00 UTC, disappeared at 16:00, which was consistent with the time of solar radiation. Figures 19, 20 and 21 present the electron density slices at different latitudes, longitudes, and altitudes in Region II at 00:00 UTC by CS. Enhanced electron density between 300 and 500 km can be seen in these three figures, and the electron density peaked at 45° latitude and at 380 km altitude. As shown in Fig. 21, two high-electron density bands in the slice at 400 km and 45° latitude were distributed on each side of the 5° longitude line. Figures 22, 23 and 24 show the time variation of the electron density slices at the 47° latitude line, the 2° longitude line, and 300 km altitude in Region II by CS, respectively. The electron density reached the minimum between 2:00 and 4:00 UTC; then, it began to increase and reach its maximum at 18:00 UTC; after that, it decreased again. The peak height decreased to a minimum at 8:00 UTC and then increased with time. In addition, the graphs show that the horizontal gradient in the longitude direction was smaller than that in the latitude direction.
In addition, CT by CS shows high computational efficiency. The computational efficiency is affected by the iterations of the algorithm and the observed quantity. In our case, the number of iterations was mostly less than 10. The observed quantity used in one calculation was about 2000, 3 epochs for Region I, 21 epochs for Region II. The algorithm ran on an ordinary desktop computer; Table 5 shows the runtime environment. Figure 25 presents the elapsed time of each reconstruction by CS; most of them were less than 10 s, and the mean values for Regions I and II were 6.1 s and 5.7 s, respectively. It should be noted that the STEC g and the dictionary were prepared in advance. Saito et al. (2017) showed that their method can produce three-dimensional electron density distributions over Japan every 15 min with a latency of about 6 min, the cost of DCB estimation was within 1 s.

Summary
Ionospheric tomography is an ill-posed problem. CS is able to reconstruct the signal through a small number of observations under certain conditions. In this study, CS was first used for ionospheric tomography. Specifically, a dictionary was constructed based on the IRI model, and the sparse representation of the ionospheric electron density was studied. Moreover, using weight matrix optimization to improve the performance of the observation matrix and the scale factor obtained through the NeQuick 2 model to retain more observations ensured the performance of the observation matrix. GPS data in two regions (I: China, II: Europe) were used to reconstruct the electron density based on this proposed method, and the results were compared with ionosonde observations. The main results of this study are summarized as follows: 1. In Region I, the peak electron density reconstructed by CS is in agreement with the ionosonde observation; the peak height is generally higher than the observation; the reason may be that the inaccuracy of IRI in Region I affects the suitability of the dictionary. In Region II, the peak electron density and height are consistent with the ionosonde observation. 2. Distribution and variation of the electron density reconstructed by CS in both regions were consistent with the actual situation. The results showed that the horizontal gradient of the electron density in the latitude direction was larger than that in the longitude direction. In addition, the electron density exhibited the local time dependence. 3. Ionospheric tomography based on CS greatly improved the resolution and computational efficiency. The average time for each reconstruction was about 6 s. For Region I, each reconstruction only uses data within 1 min. Although the calculation is conducted once every 5 min, to match the temporal resolution of the ionosonde data, it can be flexibly adjusted, such as one reconstruction once per epoch.

Data Availability
The GPS data of Region I are currently not publicly available but can be provided upon reasonable request. The data of Region II can be downloaded from http:// www. gage. upc. edu/ produ cts,

Fig. 24
Electron density sections along the 300 km height line at different times in Region II and ionosonde data can be downloaded from https:// lgdc. uml. edu/ common/ DIDBF astSt ation List.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.