Medium-scale traveling ionospheric disturbances by three-dimensional ionospheric GPS tomography

In this study, we develop a three-dimensional ionospheric tomography with the ground-based global position system (GPS) total electron content observations. Because of the geometric limitation of GPS observation path, it is difficult to solve the ill-posed inverse problem for the ionospheric electron density. Different from methods given by pervious studies, we consider an algorithm combining the least-square method with a constraint condition, in which the gradient of electron density tends to be smooth in the horizontal direction and steep in the vicinity of the ionospheric F2 peak. This algorithm is designed to be independent of any ionospheric or plasmaspheric electron density models as the initial condition. An observation system simulation experiment method is applied to evaluate the performance of the GPS ionospheric tomography in detecting ionospheric electron density perturbation at the scale size of around 200 km in wavelength, such as the medium-scale traveling ionospheric disturbances.


Background
Since late 1980s, the satellite radio tomography method has been employed to study the ionospheric electron density structure. Austen et al. (1986Austen et al. ( , 1988 first proposed the ionospheric tomography technique and measured the ionospheric total electron content (TEC) along the line of sight (LOS) from naval navigational satellite system (NNSS) to the ground-based receivers. They further successfully reconstructed the two-dimensional image of the ionospheric electron density from the onedimensional TEC data by using the technique of computerized ionospheric tomography (CIT). After that, many papers were published regarding the improvement and applications of CIT (Austen et al. 1988;Raymund et al. 1990;Yeh and Raymund 1991;Andreeva et al. 1992;Fremouw et al. 1992;Pryse and Kersley 1992;Fehmers 1994;Kunitake et al. 1995), which show the CIT has become a usable radio technique for the ionospheric electron density reconstruction (Raymund et al. 1994). Recent studies (Pryse 2003;Kersley 2005) presented that combination of the CIT with ionospheric models could help further explore the variations of global/regional ionospheric electron density structures in more detail.
Previous CIT studies usually used the high-rate NNSS observation data to reconstruct the ionospheric electron density structure (Raymund et al. 1990; Kunitake et al. 1995). This kind of low earth orbit (LEO) satellites normally has relatively high ground-track velocity. Therefore, the ionospheric electron density structure can be seen as a "snapshot" for the density inversion during the satellite passing by (within 15-30 min). In the recent decades, scientists use the ground-based global position system (GPS) receivers to derive the TEC along GPS satellites to the receivers. Because GPS-TEC data have fine spatial and temporal resolution, they are used to present the two-dimensional TEC map in the longitude-latitude plane (Wilson et al. 1992;Sardon et al. 1994;Liu et al. 1996). If we employ the CIT technique to the GPS-TEC observation data, we could obtain the altitude information of ionospheric electron density profile additionally. However, due to the slow-moving ground tracks and high elevation angle of GPS satellite, the spatial distribution of GPS-TEC observation data is not enough to reconstruct the three-dimensional electron density structure if the ground-based receivers are sparse. A dense GPS receiver network, such as Japan GPS Earth Observation Network (GEONET), is required to overcome the limitation.
Although there are regions of dense ground-based GPS receivers that could be chosen to perform GPS tomography, it is still difficult to straightforwardly solve the illposed inverse problem. This is due to the fact that the high orbit of GPS satellites (around 20,200 km) results in fewer low elevation LOS ray paths and thus less altitudinal information. The other problem is the lack of observations at the higher altitude region resulting in numerous areas that have no LOS ray passing through. Therefore, scientists employ an iterative method with an initial guess value from the ionospheric model, such as the international reference ionosphere (IRI) model, to approach the true solution of ionospheric electron density (Raymund 1995;Bilitza and Reinisch 2008;Ssessanga et al. 2015). As a result, the solutions might be sensitive to the initial guess. Therefore, it might not be suitable for the ionospheric disturbance period when ionospheric electron density may become irregular, such as the periods of solar flare or magnetic storm. In this study, an ionospheric GPS tomography algorithm with the constrained least-square method is developed and used to reconstruct the threedimensional electron density structure without using any initial model guess. An observation system simulation experiment (OSSE) is further employed to evaluate the performance of GPS tomography method during quiet time and disturbed ionospheric period.

Methods
The TEC value along the LOS from the GPS satellite to the ground-based receiver can be calculated by the following equation: in which N e = electron density, ds = length along the radio ray path, sat = location of GPS satellite, and rec = location of ground-based GPS receiver. In order to inverse the three-dimensional electron density structure from the two-dimensional TEC observation, we divide the three-dimensional space into small grids with the altitude from 80 to 20,200 km ( Fig. 1). Furthermore, we assume the electron density in each small grid has the same value. Therefore, the GPS-TEC observation along the ray path from a GPS satellite to a ground-based GPS receiver can be presented as: where A ij = length of path i in grid j, x j = electron density in grid j, and b i = TEC value along one GPS observational However, due to lack of the horizontal GPS-TEC observation paths and the numerous grids in A matrix that have no LOS ray passing through, the electron density (x matrix) cannot be calculated by straightforward solution of the least-square method. Here, we propose a modified methodology by applying a constrained condition in the least-square fitting to solve this ill-posed inverse problem.
The cost function and the solution (electron density, x matrix in 2) of constrained least-square method are showed as: in which J(x) is the cost function, which shows the normal least-square fitting at the first term of right-hand side and the constrained condition at the second term of right-hand side. W is the constraint matrix, which is a zero matrix but the value is 1 at the six neighbor grids of x. Therefore, Wx presents the total electron density differences of each grid, x j , from its six neighbor grids, x jk . C is the matrix of constraint parameter and indicates the weighting of constrain condition which depends on the altitude. Since the electron density around the ionospheric E and F regions has diurnal variations, a relatively small constraint parameter is applied from 180 to 650 km. Table 1 presents the values of the constraint parameters used in this paper.
is the regularization parameter, which is suggested by balancing the least-square term Schematic map of grid division between GPS satellite and ground-based GPS receiver from 80 to 20,200 km altitude (Chen 2012) and the constraint condition in (3) (Seemala et al. 2014). The ionospheric GPS tomography algorithm employed in this paper is similar to that described in Seemala et al. (2014), but using the fixed constraint parameters to totally remove the model dependence and further using the smaller spatial resolution in latitude and longitude to evaluate its performance in detecting relatively smallscaled ionospheric density structure.

Results for quiet ionosphere
In this paper, we choose 1241 ground-based GPS receivers from GEONET and International GNSS Service (IGS) around Japanese region for evaluating the ionospheric GPS tomography. Figure 2 indicates the location of GPS receivers and the horizontal range of tomography. The horizontal resolution is 1° in latitude and longitude within the blue area (6°N-70°N in latitude and 100°E-165°E in longitude). The black area (6°N-54°N in latitude and 120°E-154°E in longitude) has 2° horizontal resolution. Outside the black area, the resolution becomes 5°. The vertical range of tomography is set to the altitude from 80 to 20,000 km. The vertical resolutions are 20 km within 80-500 km range, 50 km within 550-900 km range, 100 km within 1000-2000 km range, and 5000 km within 5000-20,000 km range, respectively.
In order to generate the synthetic TEC observations for the OSSE, we use the ionospheric empirical model, NeQuick (Radicella and Leitinger 2001), as the model truth. Based on the LOS geometry between GPS satellites and ground-based GPS receivers, the synthetic TEC data are calculated by integrating the NeQuick electron density in each grid box along the GPS observation ray path. The synthetic GPS-TEC observations in the OSSE correspond to realistic GPS-to-receiver geometries that are collected at 10:00 UT on May 23, 2012. Furthermore, the synthetic TEC values are assumed to have no observational error in this study. Figure 3 shows the two-dimensional electron density results of GPS tomography on altitude-latitude plane at 136°E longitude (Fig. 3a-c) and latitude-longitude plane at 300 km altitude ( Fig. 3d-f ), respectively. Comparing with the OSSE model truth given by NeQuick (Fig. 3a) and tomography (Fig. 3b) electron density in the altitude-latitude plane, it can be seen that the tomography results show a similar electron density distribution with the model truth, high (low) electron density at low (high) latitude region. The reconstructed electron density error shown in Fig. 3c indicates a well reproduction of OSSE model truth. It shows that the reconstructed error is located at the range of −3.0 × 10 11 to 2.0 × 10 11 el/m 3 . Figure 3d, e shows the horizontal electron density structures at 300 km altitude on the latitude-longitude plane for OSSE model truth given by NeQuick and tomography result, respectively. It is clearly seen that the tomography results resemble the electron density distributions given by the model truth, indicating a good tomographic reconstruction. Their difference shown in Fig. 3f further suggests that the GPS tomography method can well reproduce the model truth, except for the region of 25°N-30°N and 140°E-145°E. This discrepancy might be caused by the lower LOS ray paths around this region.
The root-mean-square error (RMSE) is further used as the accuracy measurement of the tomography reconstruction. The RMSE is defined as the root-mean-square difference between the tomography electron density and the electron density of OSSE model truth, and can be described as in which N is the total number of grid boxes with the radio ray path passed, x j,tomo is the tomography electron density in each grid box, and x j,model is the model truth electron density given by NeQuick. In this case, the RMSE is around 7.97 × 10 10 el/m 3 and the average model truth electron density at 300 km altitude is around 9.83 × 10 11 el/m 3 . The RMSE value of the tomography is smaller than that of the average 300 km electron density in the ionosphere with the ratio of 8.1 %, which confirms the reliability of the inversion reconstruction using the GPS tomography algorithm proposed in this paper.
A histogram shown in Fig. 4a further presents the overall reconstructed error distribution of tomography electron densities difference with the model truth densities at each grid box in the defined tomography range. Results show that the error distributes like a Gaussian function. Around 66 % of the tomography results are within ±3.75 × 10 10 el/m 3 error range, which is ±3.8 % of the average electron density at 300 km altitude. Furthermore, we compare the tomography-reconstructed vertical density profile with the OSSE model truth at three locations. These three locations are the southwest of Japan (26°N, 128°E), the center of Japan (36°N, 136°E), and the northeast of Japan (40°N, 140°E), respectively (red points in Fig. 3f ). Figure 4b-d indicates that the reconstructed hmF2 and NmF2 are in good agreement with the model truth, especially above the hmF2. However, the reconstructed electron densities are larger than the model truth at the bottom side of F2 peak. Furthermore, the reconstructed errors are small for the NmF2 with larger error at the lower latitude region (Fig. 4b) than at the higher latitude region (Fig. 4d). The reconstructed hmF2 is located at 320 km altitude, but the model truth of hmF2 is 300 km altitude, which has ~6.67 % enhancement of 300 km altitude. This different hmF2 altitude might be caused by the vertical resolution (20-km resolution around F peak) of tomography in the study. The result, however, validates the performance of the GPS tomography algorithm to reconstruct the three-dimensional electron density from the two-dimensional TEC data without the ionospheric initial guess.

Results for ionosphere with MSTID
The traveling ionospheric disturbances (TIDs) are wavelike electron density disturbances that propagate through the ionosphere and cause wave-like TEC disturbances (Davies 1990; Kelley 2011) due to external energy input to the ionosphere and/or plasma instabilities. In midlatitude, medium-scale TID (MSTID) appears frequently. Its two-dimensional structure in airglow images and TEC disturbances were further observed by ground-based all-sky imager (Ogawa et al. 2002;Shiokawa et al. 2003;Otsuka et al. 2009) and dense GPS receiver network over Japan (Saito et al. 1998;Tsugawa et al. 2004;Otsuka et al. 2011), respectively. However, the TID-triggered threedimensional electron density structures are hardly known Fig. 2 Horizontal range of GPS tomography in this study. The latitude range is from 6°N to 70°N and the longitude is from 100°E to 165°E. The horizontal resolutions are 1° within the blue area (30°N-40°N in latitude and 130°N-140°N in longitude), 2° within the black area (6°N-54°N in latitude and 120°N-154°N in longitude), and 5° outside the black area, respectively. Red point indicates the location of ground-based GPS receiver by using 630-nm all-sky images and GPS-TEC observations. The ionospheric tomography is one of the methods to reconstruct the three-dimensional density structure. In order to evaluate whether the GPS tomography algorithm in this paper could be employed to reconstruct the ionospheric perturbations, we attempt to reproduce the three-dimensional density structure of MSTID from the two-dimensional GPS-TEC observations over Japan by performing an OSSE. Since the MSTID is defined as a medium size of TID with the horizontal wavelengths of 100-250 km (Hunsucker 1982), a density perturbation caused by electric filed during the nighttime (c.f. Otsuka et al. 2013) is applied to the quiet-ionosphere electron density for synthetic MSTID electron density as: where v′ (=25 m/s) is the drift velocity caused by the polarization electric filed. ω is the angular frequency (=2π/T, where T is the propagation period of 40 min). I and D indicate the inclination angle (=45°) of the magnetic field line and the azimuth angle (=225°) between the MSTID propagation direction and the magnetic north, respectively. The horizontal wavelength is set as 200 km to match the wavelength definition of MSTID. As a result, the upward and downward motions of electron densities are observed above 30°N latitude (Fig. 5a). Furthermore, the wavefront of electron density perturbation is in northwest-southeast (NW-SE) direction, a typical MSTID waveform alignment, as shown in Fig. 5d. Figure 5b illustrates the reconstructed electron density distribution by GPS tomography as a function of altitude and latitude at the 136°E longitude plane. The tomography result shows several clear density perturbation regions above 30°N latitude and around 300 km altitude, which is similar to the MSTID structure from the OSSE model truth (Fig. 5a). The distribution of reconstructed electron density error also has the similar pattern with the quiet time error (Fig. 3c) and similar error range of −2.0 × 10 11 to 3.0 × 10 11 el/m 3 (Fig. 5c). Figure 5d, e shows a comparison of the horizontal electron density structure by model truth and tomography reconstruction at 300 km altitude. Meanwhile, the typical NW-SE alignment of MSTIDs can be seen by the tomography reconstruction around 130-140°E longitude (Fig. 5e), similar to the perturbed density structure for synthetic MSTIDs (Fig. 5d). This comparison indicates that the proposed GPS tomography algorithm is capable of reconstructing ionospheric perturbation of 200 km scale size based on the OSSE. The RMSE is also calculated for the MSTID case. The RMSE value is around 8.61 × 10 10 el/m 3 , which is around 8.8 % of the average electron density at 300 km altitude. Comparing with the RMSE of quiet time case and MSTID case, the RMSE of MSTID case is a little larger than that of quiet time case. This might be caused by the disturbance density of MSTID. However, the tomography results still reproduce the synthetic MSTID density structure well. The distribution of overall tomography inversion error is further shown in Fig. 6a. The Gaussian distribution of error shows that around 62 % of tomography results locate within the error range of ±3.75 × 10 10 el/m 3 . The vertical density profiles at different latitudes are further presented in Fig. 6b-d. It can be seen that hmF2 and NmF2 determined by the GPS tomography algorithm are observed at almost the similar altitudes and densities as the OSSE model truth, except the NmF2 at the location of 40°N latitude and 140°E longitude (Fig. 6d). The reconstructed hmF2 at the low and high latitudes (Fig. 6b, d) is overestimated around 20 km in altitude between tomography and model truth, which might be due to the vertical resolution at the F peak.

Discussion and conclusions
Because of the lack of horizontal GPS-TEC observation paths and incomplete path length matrix (aforementioned A matrix), the previous ionospheric tomography  (Austen et al. 1986) and multiplicative algebraic reconstruction technique (MART) (Raymund et al. 1990;Ssessanga et al. 2015), usually use the iteration method with an initial guess value to reconstruct the structure of the ionospheric electron density. These iteration methods are a mathematical procedure that generates a sequence of improving approximate solutions for the inverse problems. However, the solution is sensitive to the initial guess value and becomes the same as the initial value for grids without any radio ray path passing by.
In order to resolve the limitations of the GPS tomography, the algorithm of constrained least-square method is developed and the dense ground-based GPS receiver network around Japan region on May 25, 2012, is employed to evaluate its performance. Furthermore, the reliability of this GPS tomography algorithm is validated through the OSSE during the quiet time of ionosphere. This result indicates that an appropriate three-dimensional electron density structure can be reconstructed by using this GPS tomography algorithm with the ground-based GPS-TEC observations. Then, the electron density perturbation structures with 200 km wavefront above 30°N latitudes are further applied to the model truth given by NeQuick for synthetic MSTID structure. The GPS tomography is successfully employed to reconstruct electron density distribution during the disturbed ionosphere. The tomography results have shown that the density perturbation structure can be reconstructed well. Otsuka et al. (2013) simulated the density structure of MSTID and found the different morphologies for the daytime and nighttime MSTIDs. Although the mechanism of MSTID is still not yet clear, the daytime MSTID might be caused by a gravity wave, while the nighttime MSTID might be generated by the electrodynamical coupling between the Es layer and F region (Shiokawa et al. 2003;Saito et al. 2007;Otsuka et al. 2007;Seker et al. 2008;Yokoyama et al. 2009). If the MSTID is produced by the electrodynamical coupling effect, a field-aligned structure of density should be observed because the electric field prefers to transport along the magnetic field. In this paper, we simulated the electron density structure of nighttime MSTID by Otsuka et al. (2013) and successfully reconstructed the MSTID density structure by the GPS tomography. It indicates the capability of GPS tomography for the studies of disturbed ionosphere, such as the MSTID. In the future, the developed algorithm of GPS tomography in this paper will be used to reconstruct the three-dimensional electron density structure by the real GPS-TEC observations from GEONET and further provide the near-real-time ionospheric density structure over Japan. In the present version of GPS tomography process, it needs around 30 min (including the GPS-TEC calculation and the reconstruction times) to reproduce the 3D structure of electron density. A Python version of GPS tomography is ongoing to develop to speed up the whole processes. The first result shows the time can be reduced to around 5 min, which could make the "near real time" possible after the program being completed.

Appendix: Solution of constrained least-square method
In order to solve the electron density (x matrix) from the cost function J(x) (Eq. 3), at the first step, we expand Eq. (3) as follows: where the W matrix in Eq. (7) is given as: 6C 1,1 −C 1,2 −C 1,3 −C 1,4 −C 1,5 −C 1,6 6C 2,1 −C 2,2 −C 2,3 −C 2,4 −C 2,5 −C ,6 0 . . . in which N is the total number of observation path and C jk matrix is the constrain parameter depending on the altitude, which is given in Table 1. At the next step, we assume the partial difference equation of ∂J (x)/∂x equals 0 to get the optimization solution. Therefore, Eq. (7) can be written as: Then, we move the second term of Eq. (9) to the right side of equation and can calculate the solution (x matrix) as follows: Finally, we can get the result by Eq. (11), which is the same as Eq. (4) in the text.