Introduction

The ionosphere is the upper atmosphere of the earth, where a small but significant number of the neutral particles are ionized, resulting in free electrons and ions, namely a plasma (Bust and Mitchell 2008). As a dispersive medium, the ionosphere influences the direction and the velocity of a radio wave. The effects are in proportion to the frequency and the ionospheric electron density (IED). In order to specify the IED distribution, some ionospheric theoretical and empirical models were put forward based on ionospheric observations (Anderson et al. 1998). The international reference ionosphere (IRI) is the most widely used empirical model and specifies the electron density below 2000 km (Swamy et al. 2013). Nevertheless, the IRI is kept improving especially in the topside ionosphere (Mosert et al. 2007). However, conventional measurement techniques are limited. For instance, an ionosonde can only observe the bottom side ionosphere and an incoherent scatter radar are too expensive to be set up widely. The advent of the global navigation satellite system (GNSS) enables the ionospheric total electron content (TEC) to be observed flexibly and efficiently based on group delays and phase advances of the signals (Wilson et al. 1992). With a GNSS multi-station network, computerized ionospheric tomography (CIT) can be carried out (Austen et al. 1988), where the ionosphere is divided into three-dimensional voxels traversed by many satellite-receiver rays. By effective solution method, the electron densities of the illuminated voxels can be reconstructed.

Although the IED distribution can be reconstructed by CIT from GNSS observations, several problems with this algorithm should be solved. Because some voxels are not traversed by a satellite-receiver ray, the unique solution of the electron density of the unilluminated voxel does not exist, and hence, CIT is subject to an ill-posed problem. In addition, because of the lack of horizontal rays, the non-uniform density of the rays and the limited amount of observation stations in addition to the limited receiving angles, the condition number in CIT is large, and the results are susceptible to observation errors. Therefore, CIT is affected by an ill-conditioned problem (Wen et al. 2007). In order to solve the ill-posed and ill-conditioned problems, various effective solutions were proposed. Yao et al. (2014) classified them into non-iterative reconstruction methods and iterative reconstruction methods. The most classic non-iterative reconstruction methods are constrained least square method and Tikhonov regularization, which utilize constraints to overcome the deficiency of known condition (Chen et al. 2016; Wang et al. 2016). The constraints can come from smoothing conditions, empirical models, or additional measurements. The precision of the results depends on the constraints and constraint factors. It should be pointed out that the initial model-based constraints may reduce the accuracy of the result because of the deviation between the initial value and reality. Some other scholars studied the data assimilation based on multi-source measurements or empirical models (Minkwitz et al. 2016). Other non-iterative reconstruction methods are mainly mathematical algorithms, for instance, singular value decomposition (SVD), Bayesian estimation, wavelet and neural network (Erturk et al. 2009; Norberg et al. 2015; Ghaffari Razin and Voosoghi 2017; Ma et al. 2005). Among the iterative reconstruction methods, the algebraic reconstruction technique (ART) is the most classic (Austen et al. 1988). ART produces iterative results with the initial values from empirical models or other measurements. The improved iterative reconstruction method, named multiplicative ART (MART), can obtain nonnegative results by multiplicative iteration (Raymund et al. 1990). In addition, there are improved iterative reconstruction methods such as simultaneous ART (SART) (Nandi and Bandyopadhyay 2015) and constrained ART (CART) (Wen et al. 2010).

Since the electron densities of the unilluminated voxels are not reconstructed directly, some scholars reconstructed them by smoothing and fitting (Wen et al. 2012; Yao et al. 2013). However, these methods may result in unreasonable values in the marginal voxels and areas with a large number of unilluminated voxels. Also, limited by the quantity of observable rays and their elevation angles, the coupling of the rays is weak and the result is susceptible to the observation noise. For such cases, we present a CIT algorithm based on virtual reference stations (VRSs), called VRS–CIT, to decrease the unilluminated voxels and improve the precision of the derived IED by using interpolated virtual slant TEC (STEC). This idea originated from network real-time kinematic (NRTK), which is a differential GNSS positioning technology. In NRTK, the virtual pseudorange and carrier phase observations of a VRS are generated with the observations of several nearby stations selected according to the distances and the azimuth angles (Zou et al. 2013). In geomagnetically quiet days, the precision of the interpolated ionospheric delay is better than 10 cm, which is equivalent to about 0.951 TEC unit (TECU) (Li et al. 2012). The interpolation accuracy of the double differenced ionospheric delay is even better than 1 cm in mid-latitudes and 2 cm in low latitudes (Tang et al. 2016). Wielgosz et al. (2003) demonstrated that the interpolated virtual observations could be utilized for regional ionospheric mapping. The main errors of the virtual observations come from the small-scale variation of the ionosphere and random noises of the observations (Cui et al. 2018). They can be restrained by the cutoff distance in interpolation station selection and the weight of the observations in solution, respectively. In some cases, especially in winter or autumn, medium-scale traveling ionospheric disturbance (MSTID) can cause an unignorable error (Hernández-Pajares et al. 2012). The impact of MSTID can be modeled and mitigated (Hernández-Pajares et al. 2017). Therefore, the VRS can be designed to improve the results of CIT. The observation-based constraints can improve both the solvability and the accuracy of the reconstructed IED.

The theory of VRS–CIT is introduced with emphasis on selecting the interpolation stations and the generation of the virtual observations for the VRS. After that, the performance of VRS–CIT is examined by using the data from 127 GNSS stations in 2014. An intensive study is conducted to evaluate the accuracy of the virtual observations, the decrease of the unilluminated voxels, and the reconstructed IED. The results are then summarized.

Algorithms

The influence of the ionosphere on GNSS signal can be specified by a combination of the observations from different frequencies. Then, with the ionosphere divided into several three-dimensional voxels, the IED distribution can be derived. In VRS–CIT, the interpolation stations of the VRS are selected according to azimuths and distances, and the virtual observations are interpolated from the observations of the interpolation stations.

CIT principle

A GNSS satellite broadcasts information on two spread spectrum L-band radio signals. With a GNSS dual-frequency receiver, the pseudorange and the carrier phase observations can be received:

$$P_{i} = \rho + c\left( {\delta_{t} - \delta_{s} } \right) - I_{i} + T + \varepsilon$$
(1)
$$\frac{c}{{f_{i} }}L_{i} = \rho + c\left( {\delta_{t} - \delta_{s} } \right) + I_{i} + T + \frac{c}{{f_{i} }}N_{i} + \frac{c}{{f_{i} }}w + \varepsilon^{\prime}$$
(2)

where \(i\) is 1 or 2 corresponding to \(f_{1}\) or \(f_{2}\) frequency with \(f_{1}\) = 1.57542 GHz and \(f_{2}\) = 1.2276 GHz, \(P_{i}\) and \(L_{i}\) stand for pseudorange and carrier phase observations, \(\rho\) is the actual distance between receiver and satellite, \(c\) stands for the speed of light, \(\delta_{t}\) is the receiver clock error and \(\delta_{s}\) is the satellite clock error, and they contain unequal instrumental delays in different frequency observations. \(I_{i}\) is the ionospheric correction corresponding to \(f_{i}\), \(T\) is the tropospheric delay, \(N_{i}\) is the integer ambiguity of \(L_{i}\), \(w\) is the wind-up term, and \(\varepsilon\) and \(\varepsilon^{\prime}\) are the other errors and the noises in the observations.

From the difference of the observations with different frequencies, the slant TEC (STEC) along the propagation path can be extracted:

$${\text{STEC}}_{P} = \frac{{\left( {f_{1} f_{2} } \right)^{2} }}{{40.3\left( {f_{1}^{2} - f_{2}^{2} } \right)}}\left( {P_{2} - P_{1} } \right)$$
(3)
$${\text{STEC}}_{L} = \frac{{\left( {f_{1} f_{2} } \right)^{2} }}{{40.3\left( {f_{1}^{2} - f_{2}^{2} } \right)}}\left( {\frac{{cL_{1} }}{{f_{1} }} - \frac{{cL_{2} }}{{f_{2} }} + \frac{{cN_{1} }}{{f_{1} }} - \frac{{2N_{2} }}{{f_{2} }}} \right)$$
(4)

Because the pseudoranges contain large observation noises and the carrier phases include unknown integer ambiguities, it is better to extract the STEC with high precision through carrier phase smoothing pseudorange (Jin et al. 2012):

$${\text{STEC}}\left( t \right) = \omega_{t} {\text{STEC}}_{P} \left( t \right) + \left( {1 - \omega_{t} } \right)\left[ {{\text{STEC}}\left( {t - 1} \right) + \,{\text{STEC}}_{L} \left( t \right) - {\text{STEC}}_{L} \left( {t - 1} \right)} \right] (t > 1)$$
(5)

where \(t\) stands for the epoch, and \(\omega_{t}\) is the weight depending on the elevation angle. \({\text{STEC}}\) is equal to \({\text{STEC}}_{P}\) when \(t\) is equal to 1. Based on the ionospheric single-layer model that assumes all the ionospheric electrons concentrate on an infinite thin single layer, the instrumental biases are determined using the self-calibration of the observation error process (Bishop et al. 1996). In this method, a conjunction between two satellites is observed when they arrive at a certain area in the height of the ionospheric single layer at the same moment. For the conjunction, the same vertical TEC is assumed to be measured from each satellite. In the study, the configuration used for the definition of near conjunction is 1° in latitude and longitude.

In CIT, a day is divided into a series of time periods and a time period is usually 30 min. The ionosphere is divided into several three-dimensional voxels. The electron density value in a voxel is assumed to be identical during a time period and potentially changes over different time periods. The STEC is regarded as the line integral of the IED along the propagation path \(i\) and is represented as a finite sum over the voxels \(j\):

$${\text{STEC}}_{i} = \mathop \sum \limits_{j = 1}^{n} A_{ij} {\text{IED}}_{j} + \varepsilon_{i} \left( {i = 1,2, \ldots ,m} \right)$$
(6)

where \(A\) is the intercept matrix of the rays in the voxels and \(\varepsilon_{i}\) is an error that includes the instrumental biases and other discretization errors. Using all the observations, the matrix equation in CIT can be constructed:

$${\text{STEC}}_{m \times 1} = A_{m \times n} {\text{IED}}_{n \times 1} + \varepsilon_{m \times 1}$$
(7)

As explained above, Eq. (7) is ill-posed and ill-conditioned. It can be solved with the algorithm of MART (Raymund et al. 1990):

$${\text{IED}}_{j}^{{\left( {k + 1} \right)}} = {\text{IED}}_{j}^{\left( k \right)} \left( {\frac{{{\text{STEC}}_{i} }}{\langle{{\text{IED}}^{\left( k \right)} ,a_{i} }\rangle}} \right)^{{{\raise0.7ex\hbox{${\lambda_{k} a_{ij} }$} \!\mathord{\left/ {\vphantom {{\lambda_{k} a_{ij} } {a_{i} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${\|a_{i}\| }$}}}}$$
(8)

where \(a_{i}\) stands for the \(i\)th row vector in matrix \(A\) and \(\lambda_{k}\) is the relaxation factor with 0 < \(\lambda_{k}\) < 1. After the iteration converging, the IED distribution can be reconstructed with minimum residuals.

Interpolation station selection

In VRS–CIT, latticed VRSs are designed to improve the distributions of the satellite-receiver rays. For each VRS, several interpolation stations are selected from the nearby observation stations to generate virtual observations. To restrain the noise, the selection is on the basis of the distance. Meanwhile, in order to improve the geometric distribution, the azimuth angles are taken into consideration.

Figure 1 presents a schematic of VRS in the center of a circle, with three and six observation stations included in the circle in the left and right panels, respectively. The radius of the circle defines the cutoff distance. The determination of the cutoff distance is a compromise between available VRS and interpolation error. In this study, we set the cutoff distance as 150 km. Zou et al. (2013) revealed that the difference of STEC between stations was usually less than 1 m, and the maximum interpolation difference was only 1.4 cm for a cutoff distance of 150 km. In the left panel, three dashed rays with 2π/3 separations are plotted. The three interpolation stations are selected by turning the rays clockwise until each sector contains one observation station. In the right panel, with the selected interpolation stations and the rotation angle in the previous step recorded, three other interpolation stations are selected likewise with the initial angle being the rotation angle plus π/3. For a VRS, if no interpolation station is selected, this VRS should be deleted. Then the virtual observations are generated by the STEC of the six stations or the three stations.

Fig. 1
figure 1

Diagram of interpolation station selection for a VRS. The left panel and the right panel represent the selection of three interpolation stations and three additional interpolation stations, respectively. The center represents the VRS and the radius is the cutoff distance. The solid lines and the dashed lines are the rays before and after rotation, respectively. The hollow circles are the selected interpolation stations

Full size image

Virtual observation generation

The virtual observations of VRSs are interpolated from the observations of nearby six or three stations. For virtual observation generation with six stations, multi-quadric surface fitting (MSF) is utilized. When there are not enough stations available for MSF, triangular linear interpolation (TLI) is chosen.

Multi-quadric surface fitting

MSF is usually utilized to interpolate the data in a multi-dimensional space, and the theoretical basis is that any smooth mathematical surface can always be approximated with any precision by the sum of a series of regular mathematical surfaces. In fact, MSF is a barycenter-based interpolation and the distance is redefined. The weight of each interpolation station is determined by the distance to the barycenter of all the interpolation stations. Because the distributions of the interpolation stations are usually inhomogeneous, the barycenter-based interpolation has better stability and lower interpolation error.

In this method, multiple quadric surfaces with centers or vertexes located at the interpolation stations are generated. Then the virtual observation is computed by the linear combination of the surfaces. First, the weights of the interpolation stations are determined by STEC values and the distances between two stations from (9). For six stations, seven equations result to solve the unknowns \(c_{j}\) and \(b\). Then the virtual observation of VRS is computed by the solved weights and the distance between VRS and each interpolation station from (10) (Al-Shaery et al. 2011):

$$\left\{ {\begin{array}{*{20}l} {{\text{Obs}}_{i} = \mathop \sum \limits_{j = 1}^{n} c_{j} g_{ij} + b \left( {i = 1,2, \ldots n} \right)} \hfill \\ {\mathop \sum \limits_{j = 1}^{n} c_{j} = 0} \hfill \\ \end{array} } \right.$$
(9)
$${\text{Obs}}_{0} = \mathop \sum \limits_{j = 1}^{n} c_{j} g_{0j} + b$$
(10)
$$\left\{ {\begin{array}{*{20}l} {g_{ij} = \left[ {\left( {x_{i} - x_{j} } \right)^{2} + \left( {y_{i} - y_{j} } \right)^{2} + a^{2} } \right]^{1/2} } \hfill \\ {g_{ij} = \left( {x_{i} - x_{j} } \right)^{2} + \left( {y_{i} - y_{j} } \right)^{2} } \hfill \\ {g_{ij} = \left[ {\left( {x_{i} - x_{j} } \right)^{2} + \left( {y_{i} - y_{j} } \right)^{2} } \right]^{1/2} } \hfill \\ \end{array} } \right.$$
(11)

where \(c_{j}\) is the weight of each quadric surface and the sum of the weights is assumed to be zero. The constant term \(b\) is designed for the non-bias property, \(n\) is the number of the quadric surfaces, here \(n = 6\), and \(g_{ij}\) is the distance parameter, which can be a hyperboloid surface, a paraboloid surface, or a conic surface showed in (11). Since for the hyperboloid surface, the optimum constant term \(a\) is hard to determine and for the paraboloid surface there are computational problems, a conic surface is usually chosen.

Triangular linear interpolation

TLI can generate the virtual observations of a VRS in a triangle network with three stations located at the vertexes according to the distances. As depicted in Fig. 2, with the coordinates of interpolation stations \(A\), \(B\),\(C\) and VRS \(E\), the interpolation value can be computed as follows (Mannucci et al. 1998):

$${\text{Obs}}_{E} = \frac{{{\text{ED}}}}{{{\text{AD}}}}{\text{Obs}}_{A} + \frac{{{\text{AE}}}}{{{\text{AD}}}}\left( {\frac{{{\text{DC}}}}{{{\text{BC}}}}{\text{Obs}}_{B} + \frac{{{\text{BD}}}}{{{\text{BC}}}}{\text{Obs}}_{C} } \right)$$
(12)
Fig. 2
figure 2

Diagram of TLI. The solid circle stands for the VRS, and the hollow circles stand for the interpolation stations. The point \(D\) is the intersection of lines \(AE\) and \(BC\)

Full size image

It should be noted that the VRS out of the triangle network ought to be eliminated.

Test and validation

In the following subsections, we examine in detail the performance of VRS–CIT with dual-frequency data from the global positioning system (GPS) on April 29, 2014, which is a geomagnetically quiet day in a solar maximum year. Then, an adaptability analysis is conducted with twelve geomagnetically quiet days of 21–23 in March, June, September and December in 2014. In order to have a comprehensive assessment of VRS–CIT, CIT is conducted with only observation stations (OBS–CIT) and then with the latticed VRSs added (VRS–CIT) or with the horizontal Laplacian constraint added (CCIT) (Wen et al. 2010).

Data and processing method

The data used to examine the performance of VRS–CIT are from 127 observation stations from the GNSS earth observation network system (GEONET) and the international GNSS service (IGS). Five additional observation stations from GEONET are selected for validation. As shown in Fig. 3, the top panel illustrates the distributions of the observation stations, and the bottom panel presents VRS distribution.

Fig. 3
figure 3

Distributions of the observation stations and VRSs. The meshes with the dashed lines in the top panel represent the main area in the unsolvability analysis

Full size image

The epoch interval of the data is half a minute, and the cutoff elevation angle is set as 20° to restrict the interpolation error and the error from the multipath effect. The density of the VRS is 0.5° × 0.5°, as shown in the bottom panel. Before the tomography, the instrumental biases are eliminated. In CIT, the initial values come from IRI-2016, and MART is utilized. The reconstructed area covers 24–46° N in latitude and 122–146° E in longitude. The altitude ranges from 100 to 1000 km. The spatial resolution is 2°, 2°, and 25 km in latitude, longitude, and altitude, respectively. The time resolution is half an hour. In consideration of the interpolation error, the weights of real observations and virtual observations by MSF and TLI are set as 1, 1/2, and 1/3, respectively. An assessment proceeds to be conducted on the decrease of the unilluminated voxels and the improvement of the precision. To specify the precision, the indexes of the mean absolute error (MAE) and the root mean square error (RMSE) are computed as follows:

$${\text{MAE}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left| {y_{d} - y_{r} } \right|$$
(13)
$${\text{RMSE}} = \sqrt {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {y_{d} - y_{r} } \right)^{2} }$$
(14)

where \(N\) is the amount of the derived values, \(y_{d}\) and \(y_{r}\) are the derived and real values, respectively.

Virtual STEC assessment

Five validation stations depicted in the top panel of Fig. 3 are selected to evaluate the accuracy of the virtual observations generated by MSF and TLI. With the real observations of the validation stations taken as references, MAE and RMSE of the virtual STEC are computed from (13) and (14).

In Fig. 4, the MAE values of the virtual STEC generated by MSF and TLI fluctuate between 0.7 TECU and 1.7 TECU in the day. The average MAE values of MSF and TLI in the day are 1.0587 and 1.0779, respectively. In consideration of the observation errors in the real observations, the result accords with that of Li et al. (2012). During the day, except for 12–13 UT, the MAE values of MSF are lower than TLI. In theory, only if there are three stations located in different directions and close to VRS, the accuracy of TLI is close to that of MSF. TLI is less stable than MSF because of the inhomogeneity of the interpolation stations. However, MSF is more susceptible to ionospheric activity because further stations are used. The fluctuation of MAE and RMSE values is related to the diurnal variation of the ionospheric plasma densities. In 20–21 UT, the MAE and RMSE values significantly increased because after local sunrise, the increasing ionospheric horizontal gradient in the interpolation region causes large interpolation error. Therefore, with the error restrained by the cutoff distance, MSF and TLI can be utilized to generate the virtual observations.

Fig. 4
figure 4

MAE and the RMSE of the virtual observations for the five validation stations. The real observations of the validation stations are taken as references

Full size image

Unsolvability analysis

Limited by the quantity and inhomogeneity of the satellite-receiver rays, some voxels are not illuminated. The electron densities of the unilluminated voxels are not solved and keep equal to the initial values in the test. To investigate the decrease of the unsolvable voxels with the VRS added, relatively small spatial resolution is chosen for tomography. With the spatial resolution mentioned above, the reconstructed area yields a total of 4752 inversion variables. The spatial resolution is close to that of Yang et al. (2017); In the results, the unsolvability of the electron densities for the voxels is also close to theirs, which is about 30%.

In the top panel of Fig. 5, for OBS–CIT there are still nearly 15% of the whole voxels unsolvable around 2 UT and 13 UT. These unsolvable voxels are mainly located at the edge of the reconstructed area. The selected observation stations range from the northeast to the southwest of the reconstructed area. The voxels located at northwest and southeast may be unsolvable, especially at the bottom-side ionosphere, where the rays are relatively sparse. Then in the bottom panel, for the main area depicted in the top panel of Fig. 3, there are very few unsolvable voxels around 2 UT and 13 UT. During 6–7 UT, the unsolvable voxels are over 30% for the whole area and 10% for the main area. It is related to bad distributions of the rays during this time period.

Fig. 5
figure 5

Unsolvability of the electron densities for the voxels in CIT. The black points denote the decline rate of unsolvability in VRS–CIT comparing to that in OBS–CIT

Full size image

Comparing to OBS–CIT, the unsolvability of VRS–CIT declined by 4–12% for the whole area, and the rate is even larger for the main area. In the bottom panel, the reason for the unstable decline rate is that few voxels in the main area are unsolvable during some time periods. Disregarding the extreme cases, the most obvious decline rate of unsolvability is about 70% at 15 UT. With the distributions of the stations improved and satellite-receiver rays increased by latticed VRSs, some voxels are traversed by virtual rays and the electron densities become solvable. Therefore, the VRS is significant in the solution of ill-posed and ill-conditioned problems. Even though the majority of the improvements are located in the main area, the unreasonable smoothing and fitting in the marginal voxels and the area with a large number of unilluminated voxels can be avoided. Furthermore, if some observation stations can be added at the edge of the reconstructed area, the decline rate can rise further.

Vertical profile validation

In order to assess the accuracy of the IED derived by VRS–CIT, the vertical profiles of IED from radio occultation are taken as references. The FORMOSAT-3/COSMIC constellation contains six satellites with 30° separations and keeps an orbital period of 100 min. The orbit maintains an altitude of 800 km for the period of this study (Farzaneh and Forootan 2018). The second-level data ‘ionprf files’ describe the vertical profiles of IED. The data provided by COSMIC data analysis and archive center (CDAAC) have a reported accuracy of \(1 \times 10^{10} - 1 \times 10^{11} \;{\text{e}}/{\text{m}}^{3}\) (CDAAC 2013).

In Fig. 6, the vertical profiles of IED reconstructed by CIT are compared with the profiles from COSMIC. In the top panel, the COSMIC line is the only one shown above 450 km. The reason is that the tangent point of COSMIC is out of the reconstructed area during this time period and the IED cannot be derived by CIT. Around 250 km, the IED by OBS–CIT is closer to the IRI profile than VRS–CIT, while the IED by VRS–CIT is closer to the COSMIC profile than OBS–CIT. It demonstrates that with observation stations only, the electron densities are close to the initial values for some voxels, but with the VRS added the electron densities of these voxels become more accurate. The IED by CCIT is closer to the IRI profile, compared with that by VRS-CIT. The reason is that the constraint from the initial model is enhanced for CCIT. In the bottom panel, similar results can be derived around 300 km. Therefore, the IED by VRS–CIT is more accurate, especially around the height of IED maximum. With the COSMIC profiles taken as reference, the MAE of the CIT profile is computed. Comparing to OBS–CIT, the MAE of VRS–CIT decreases by 6.88% and 8.43% for the two IED profiles, respectively. Therefore, the precision of derived IED is improved with the VRS added.

Fig. 6
figure 6

Vertical profiles of IED reconstructed by CIT. The top panel and the bottom panel are corresponding to two respective occultation events. The initial values of electron densities offered by IRI-2016 are illustrated. The title shows the local time at the moment of IED maximum in the occultation event

Full size image

Estimated STEC validation

In theory, the IED of each voxel in the reconstructed area can be derived by CIT. Nevertheless, in fact, the electron densities of the unilluminated voxels are equal to the initial values in the test. The intercept of a ray in each voxel can be calculated by geometry. With the IED and the intercepts, the STEC of the validation stations can be estimated by the line integrals and proceed to be compared with the STEC extracted from the real observations of five additional validation stations.

In Fig. 7, the local peak values of MAE and RMSE are consistent with those of the unsolvability, for instance, at 3 UT, 7 UT, and 9 UT. It demonstrates that fewer solvable electron densities result in lower precision of the reconstructed IED distribution. The MAE and RMSE of CCIT are larger than those of OBS–CIT, especially at 6 UT. The reason is that at local noon, the initial model-based constraints cause a larger error because of the high plasma densities and large ionospheric horizontal gradient. With the VRS added, the accuracy of the reconstructed IED is improved. The maximum decline rates of MAE and RMSE are 17.24% and 33.81% at 7 UT, respectively. With some voxels traversed by the virtual rays, more values are not equal to the initial values. Because of the spatial correlation of the ionosphere, the virtual STEC is close to the realities, and the IED of the voxel traversed by virtual observation is more accurate than the initial values. In addition, the virtual observations enhance the coupling of the rays, and hence, the utilization rate of the initial values is improved. As a result, the MAE and the RMSE of VRS–CIT are lower than OBS–CIT.

Fig. 7
figure 7

MAE and the RMSE of the STEC values, reconstructed by the line integrals of the electron densities. The STEC values computed by the observations of five additional validation stations are taken as references

Full size image

However, in only a few time periods, the MAE and the RMSE rise slightly. It may be related to noise in the virtual observations. The virtual observation is generated with the linear combination of the observations and the noise may be amplified. Therefore, reasonable weighting for virtual observations is necessary.

Latitude-altitude image analysis

In the bottom-left panel of Fig. 8, the IED image of IRI illustrates that the IED increases with the decreasing latitude at mid-latitudes. But with the altitude rises, the IED first increases and then decreases with the maximum located at the height of 300–400 km. According to earlier discussion, the IRI is an empirical model and has discrepancies with the actual ionosphere. Therefore, this image represents the ionosphere in the ideal condition with no disturbance but can provide initial values to solve the ill-posed and ill-conditioned problems.

Fig. 8
figure 8

Latitude-altitude images of IED reconstructed by CIT. The longitude is 140° and the time is 3 UT. The bottom-right panel is the difference between the top-left panel and the top-right panel. The unit of IED is e/m3

Full size image

With the initial values from IRI-2016, the IED images reconstructed by CIT in the top panels are close to IRI in the bottom left panel. But at the height of about 200–500 km, the IED image reconstructed by CIT is larger than IRI in the region of 30–40° N while smaller than IRI in other latitudes. It demonstrates the discrepancies between the observational ionosphere and IRI. Comparing OBS–CIT with IRI, there is an abnormal peak of IED around 31° N and 300 km in the top-left panel, and the value is equal to that of the IRI. The reason is that there are few observation stations in the low latitude region around 140° E meridian. Therefore, for OBS–CIT the voxel corresponding to the peak is unilluminated and the IED is equal to the initial value, but the electron densities of the adjacent voxels can change with iteration. Since the electron densities of these voxels are larger than those of the adjacent voxels, the IED peak appears even it is abnormal. With the VRS added, the voxels are traversed by virtual rays and the abnormal peak disappears in the top-right panel. Therefore, the IED image reconstructed by VRS–CIT is more reasonable.

In the bottom-right panel, the difference of IED image reconstructed by OBS–CIT and VRS–CIT is relatively large in the region of 30–40° N while small in other latitudes. Because around 140°E meridian in the region of 40–45° N the observation stations are sufficiently dense, and in the region of 25–30° N the observation stations are sparse and the VRSs are very few, the effect of VRS is not apparent in these regions. In contrast, around 32°N the observation stations are relatively sparse while some VRSs are added, so the abnormal peak disappears and the IED image becomes more reasonable.

Adaptability analysis

An adaptability analysis is conducted with twelve geomagnetically quiet days of 21–23 in March, June, September, and December of 2014. The IED distributions are reconstructed by OBS–CIT, VRS–CIT, and CCIT, respectively. Then the MAE and RMSE values of virtual and estimated STEC are computed from (13) and (14).

Table 1 shows the MAE and RMSE values of the virtual STEC generated by MSF and TLI for the 12 days. The MAEs of both MSF and TLI clustered around 1 TECU, with a minimum value of 0.8354 TECU on June 23 and a maximum of 1.2011 TECU on December 23. The RMSEs of both MSF and TLI ranged between 1.1899 TECU on June 23 and 1.6810 TECU on December 23. The MAE and RMSE values in March and December are close to those of April 29, while they are smaller in June and September. In addition to larger errors on December 21–23, it is also worth noting that the MAE and RMSE values of virtual STEC of MSF are larger than those of TLI for the three days. Studies with mid-latitude stations showed that MSTID frequently occurs at local day time in winter, which increases the ionospheric horizontal gradient (Kotake et al. 2006). The larger interpolation error with 6 stations could be caused by the TEC gradients with a spatial scale-size of a few hundred kilometers due to MSTID.

Table 1 MAE and the RMSE of the virtual STEC values for the five validation stations in the 12 days in 2014
Full size table

Table 2 lists the MAE and the RMSE of the STEC values reconstructed by the line integrals of the derived electron densities for the 12 days. The MAE and RMSE values of the estimated STEC of VRS–CIT are less than those of OBS–CIT and CCIT for all 12 days. The MAE and RMSE values of CCIT can be larger than those of OBS–CIT due to the deviation between the constraint from the initial model and reality. Comparing VRS–CIT in Table 2 with MSF in Table 1, there is a positive correlation between VRS–CIT results and interpolation error. The maximum differences in MAE and RMSE are about 0.5610 TECU and 0.7683 TECU, respectively.

Table 2 MAE and the RMSE of the STEC values, reconstructed by the line integrals of the electron densities in the 12 days in 2014
Full size table

Summary

We presented a CIT algorithm based on VRSs to decrease the number of unilluminated voxels and improve the precision of the reconstructed IED distribution. In VRS–CIT, latticed VRSs are added at the nodes of grids with a 0.5° × 0.5° resolution in longitude and latitude. The virtual observations are generated by MSF or TLI with the observations of six or three nearby stations selected according to the azimuth angles and the distances. The IED distribution is reconstructed by MART and then compared to OBS–CIT and CCIT using the GPS dual-frequency data from GEONET and IGS on April 29, and another 12 days distributed in 4 seasons in 2014. With a detailed assessment on the unsolvability and the precision of derived IED, the results show that:

  1. (1)

    The precision of the virtual observations is about 1 TECU. With more voxels traversed by the virtual rays, compared to OBS–CIT, the unsolvability of VRS–CIT declined by 4–12% for the whole region, and for the main area the improvement could be up to 70%.

  2. (2)

    Taking two IED profiles from radio occultation as reference measurements, the IED profile of OBS–CIT was closer to the initial values obtained from IRI-2016 than VRS–CIT, and the IED profile of VRS–CIT was closer to the references than to OBS–CIT. With the VRS added, the MAE decreased by 6.88% and 8.43% for the two IED profiles, respectively. The most obvious improvements occurred around the height of IED maximum.

    In VRS–CIT, the decrease in the number of unilluminated voxels and the improvement of the satellite-receiver ray coupling made some derived IED values closer to the realities than the initial values. With the observations from five additional validation stations taken as references, the MAE and the RMSE of the estimated STEC corresponding to VRS–CIT are reduced. The maximum decline rates were 17.24% and 33.81%, respectively.

    For the latitude-altitude images of IED, an abnormal peak caused by unilluminated voxels at low latitude along 140°E meridian disappeared with the VRS added.

  3. (3)

    In an adaptability analysis with 12 days, the interpolation errors of MSF and TLI are about 1 TECU and less affected by ionospheric activity. With the VRSs added, the precision of CIT is improved.

In summary, with the VRS added, more IED values of the voxels in the reconstructed area can be solved, and the coupling of the satellite-receiver rays is improved, and hence, the IED distribution can be reconstructed more reasonably.