A New Ionosphere Tomography Algorithm With Two-Grid Virtual ...
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
2373
A New Ionosphere Tomography Algorithm With Two-Grid Virtual Observations Constraints and Three-Dimensional Velocity Profile Yibin Yao, Jian Kong, and Jun Tang
Abstract—Ionosphere tomography is a typical ill-posed problem, and using the ionosphere priori information as the constraints to improve the state of normal equation is an effective approach to solve this problem. In this paper, we impose priori constraints by increasing the virtual observations in n-dimensional space. Then, after the inversion region to be gridded, we can form a stable structure between the grids with loose constraints, which greatly improve the state of normal equation. Based on that, to obtain the real-time velocity information of ionosphere electron density, we introduce the grid electron density velocity parameters, which can be estimated with electron density parameters simultaneously. The authors use the new algorithm and the global navigation satellite system data in Europe to inverse the 12 ionosphere electron density and velocity images on August 15, 2003, which reflects the ionosphere changes of electron density and velocity in whole day; in addition, we compare the results with the traditional algorithm multiplicative algorithm reconstruction technique’s results. Furthermore, we compare the results with the changes time series of plasma frequency observed by ionosphere ionosonde among different layers, and many types of analysis and comparison verify the effectiveness and reliability of the new algorithm. The related research provides a new way for the real-time detection and prediction of ionosphere changes. Index Terms—Computerized ionosphere tomography, grid constraints, three-dimensional ionosphere velocity image, virtual observations.
I. I NTRODUCTION
A
S THE important part of the Earth space environment, ionosphere about 60–1000 km above the ground is the critical protective layer for human survival and development. Using different kinds of modern technologies to detect and analyze environment changes of the Sun–Earth space, particularly the Earth’s space, and then to investigate their basic structure and variation helps us not only to improve accuracy of the velocimetry, location, timing, precision communications,
Manuscript received April 17, 2014; revised July 27, 2014; accepted September 8, 2014. This work was supported in part by the National Natural Science Fund of China under Grant 41274022 and Grant 41174012 and in part by the 863 Program under Grant 2013AA122502. This paper has supplementary downloadable material available at http://ieeexplore.ieee.org, provided by the authors. This includes details about the basic principle of the new algorithm and verification of the basic idea with simulated data. This material is 214 kB in size. The authors are with the School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2014.2359762
and navigation but also to study and find out the relationship between the upper atmospheric layers. Particularly, it has important scientific significance in the study of the mechanism of global ionosphere disturbances and irregular changes. Ionosphere tomography is a typical ill-posed problem. Such problems as fewer observations in the edge of the inversion region and the observations from global navigation satellite system (GNSS) constellation are nearly vertical lead to the condition number of the normal equation in inversion algorithm being too large, and then, the estimation solution is unstable with the presence of observation noise. The commonly used reconstruction algorithms are algorithm reconstruction technique (ART), multiplicative ART (MART), and simultaneous iteration reconstruction technique (SIRT) [2], [13], [14]. The most intuitive idea to solve the ill-posed problem is increasing the observations, such as, currently, the ionosphere inversion algorithms based on multisource data [1], [5], [12], [16], [18], [19]. On the other hand, introducing the priori information of ionosphere electron density to improve the state of normal equation also can overcome the ill-posed problem. Hobiger et al. [7] used the priori constraint with SIRT to inverse ionosphere electron density and then compared the results with SIRT. After imposing the priori constraints, the iterative convergence speed and accuracy are improved. In contrast to the lack of constraint in the edge area, Liu [8] proposed an improved SIRT algorithm with priori constraints and verified its superiority. In recent years, Garcia and Crespon, Hobiger et al., Wen and Yuan, and Wen and Wang have proposed several corresponding solutions to ill-posed problems in ionosphere tomography due to insufficient projection data [6], [7], [20], [21]. However, it needs to be noted that the ill-posed problem of normal equation is still the critical problem in ionosphere tomography algorithm. Ionosphere anomaly detection, particularly the related ionosphere anomalies and ionosphere scintillation before the natural disasters (tsunamis, earthquakes), is a more popular research direction [15]. However, most of the research are confined in studying the statistical analysis of the spatial location and time position of the anomalies, which is based on the 2-D ionosphere images [3], [4], [8], [11], [17], [22]; there are fewer research about the physical mechanisms and theoretical models of the anomalies based on 3-D ionosphere images. The establishment of physical model needs the reliable and adequate ionosphere change information. The existing information for model building is insufficient to provide further information,
0196-2892 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
2374
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 1. Schematic of the computation domain of tomography.
and ionosphere electron density information in real time (ionosphere changes velocity) has an important role in the physical model establishment. With the 3-D ionosphere electron density images, the 3-D ionosphere changes velocity images can provide reliable ionosphere gradient information in each layers for some joint inversion algorithms and then improve the accuracy of ionosphere model and positioning [13], [23]. On the other hand, the 3-D ionosphere changes velocity information is conducive to ionosphere prediction and forecasting in three dimensions, which is of great significance to scientific research and engineering applications based on ionosphere changes. II. BASIC P RINCIPLES OF I ONOSPHERE E LECTRON D ENSITY R ECONSTRUCTION BASED ON THE P IXEL L EVEL The electron density reconstruction is a typical inverse problem, which is based on the observed slant total electron content (STEC) to reconstruct the 3-D distribution in formation of ionosphere electron density. As shown in Fig. 1, after the inversion region to be gridded in three dimensions, GPS satellite signal propagation paths can be taken approximately as a straight line. → Select an appropriate set of basic functions b(− r ) to model the ionosphere electron density in the inversion region, i.e., → N e(− r , t) ≈
J
→ xj (t) · bj (− r)
(1)
j=1
where J represents the number of selected basis functions, xj (j = 1, 2, . . . , J) represents the coefficients of the basis → function, N e(− r , t) is the ionosphere electron density on the → signal propagation paths, t is the observation time, and − r is the ray path vector. Then, the STEC of each ray path can be expressed as follows: STECi (t) =
J
→ xj (t)bj (− r ) ds
s j=1
=
J j=1
xj (t) s
→ bj ( − r ) ds,
i = 1, 2, . . . , M
(2)
where M represents STEC number (the GPS ray path num the → r ) ds, considering the GPS measurement ber). Let Aij = l bj (− noise, then the STEC on GPS signal propagation paths can be expressed as follows: STECi =
J
Aij xj + εi .
(3)
j=1
The formula (3) can be expressed by the vector equation as in formula (4), i.e., L = STEC = Ax · x + ε
(4)
where L is the observation vector. STEC is the column vector of STEC. Ax is the design matrix. x is the column vector composed of the coefficients of basic function. In the ionosphere tomography, x is the ionosphere electron density. ε is the column vector of the observation noise. III. N EW I ONOSPHERE T OMOGRAPHY A LGORITHM W ITH T WO -G RID C ONSTRAINTS AND T HREE -D IMENSIONAL V ELOCITY P ROFILE A. Loose Constraints Among Grids In the process of GNSS-based electron density inversion in pixel level, on one hand, due to the lack of observational data, there are no rays that pass through some grids, and on the other hand, the GNSS observations rays are almost vertical; there is limited horizontal information carried by the observation rays. Affected by the two reasons, the condition number of the normal equation matrix in parameter estimation is very large; thus, the inversion of matrix is unstable and inaccurate. There are two main ways to solve this problem. First, with the maturity and development of various ionosphere detection satellites (DORIS, occultation, satellite altimetry), the observational data types that can be used are gradually increasing. The joint inversion combining multisource data is an effective way to solve the preceding problems. In addition, although the ionosphere is affected by the Sun, geomagnetic activity, and other factors, as the important part of Sun–Earth environment, the ionosphere changes reflect the electron content change, which is the electron moving and transmission process. That
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
2375
Fig. 3. Schematic of plane constraint and altitude constraint.
Fig. 2.
Cooperation of the two different ways to add virtual observations.
means that the electron density values among different grids are correlated with each other in time and space scales. Thus, based on this, adding the constraints among the grids and then improving the state of normal equation is another effective solution. There are many ways to add constraints, such as increasing parameters constraint equations. However, the more simple and practical method is to add constraints by increasing the virtual observation to form the virtual normal equation. According to the external auxiliary information [such as global ionosphere map (GIM) and international reference ionosphere (IRI)], simulating virtual station location, virtual satellite position, and virtual observations and then combining the virtual and real observations to consist the normal equations simultaneously improves the state of the equation. Increasing the virtual observations can solve the problem of insufficient observed data in some extent, but it still cannot solve the problem that the observation data are nearly vertical, because the virtual observation values are still lacking the constraint in the horizontal direction. Thus, after adding the virtual observation, the condition number of normal equation is still large. Based on the existing algorithm, in this paper, we propose a new way to add the virtual observations, in order to increase the horizontal constraint in the normal equation and improve its state. The traditional method [see Fig. 2(a)] to add virtual observation is to simulate station location, satellite positions, and STEC; and this method is implemented in 3-D space. If we do it again in n-dimensional space, a ray from the satellite can pass through any two arbitrary grids to reach stations in
n-dimensional space [see Fig. 2(b)]. This new way to add virtual observations cannot only ensure that each grid can be passed through by virtual rays but also increase the horizontal constraint among the grids by designing the virtual rays. Then, combine the virtual observations in n dimension and the real observations to form the normal equation for parameter estimation. Compared with the traditional method, the state improvement to the normal equation of the new approach is very clear. However, if we simulate a virtual ray between every two grids, the constraints among the grids will be too tight; in case the real observations are insufficient, the grid system with tight constraints cannot be corrected rationally to keep with the real observations. Considering that the ionosphere inversion algorithms are mostly based on altitude and horizontal direction, in this paper, we add virtual observations on these two levels based on this characteristic, which is shown in Fig. 3. • horizontal constraint C: add the plane constraints in horizontal direction; • altitudinal constraint H: add the proportion constraints in altitude direction where H and C are determined by the initial value of the grid electron density. To add the H and C constraints to the estimation, first, select the grids Gi and Gj on the same altitude or at the same site orderly, and then, obtain the two grids’ electron densities N eGi and N eGj according to the IRI model. Based on this, simulate two intercepts IGi and IGj crossing the two grids artificially; here, we need to refer to the real grid size and the size range of the actual intercepts. Then, a virtual observation ¯l can be obtained with the formula (5), i.e., ¯l = IG • N eG + IGj • N eGj . i i
(5)
After forming the virtual observations among all the two grids on the same surfaces and the same sites with formula (5), ¯ i.e., we can get the virtual observations matrix L, ⎛¯ ⎞⎛ ⎞ I G1 • N e G1 + I G2 • N e G2 l1 . ⎟⎜ . ⎟ ¯=⎜ L (6) ⎝ ⎠⎝ ⎠. . . ¯ln IGn−1 • N eGn−1 + IGn • N eGn
2376
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Then, form the virtual observations to the normal equation, i.e., ⎞⎛ ⎞ ⎛ x1 0 0 IG1 IG2 . ⎟ . . . . ⎟ ⎜ ⎜ ¯=⎝ L ⎠⎝ ⎠ . . . . . xn 0 0 IGn−1 IGn ⎛ ¯ A11 ⎜ . =⎝ . 0
A¯12 . . 0
0 . . A¯nn−1
⎞⎛ ⎞ x1 0 . ⎟⎜ . ⎟ ⎠ ⎝ ⎠ = A¯x x. (7) . . A¯nn xn
The simultaneous equations can be obtained with formula (4) and formula (7), i.e., Ax L PL 0 (8) ¯ = A¯x x P = L 0 PL¯ where P is the weight matrix, PL is the weight matrix for the real observations, and PL¯ is the weight matrix for the virtual observations. In this paper, we set PL /PL¯ = 100/1. After adding these two constraints, the inversion area is formed as a stable structure with loose constraints among grids. Then, if there is a real ray passing through the area, with greater weight on the real observation, then the grid the rays pass through will selfcorrect to match real observation. Due to the plane and altitude constraints among the grids, the corrected grids will trigger the other grids in the inversion area to self-correction order to comply optimally with the real observation data; eventually, the whole inversion region achieves optimal correction based on actual observations. (More details about the basic principle of the new algorithm can be found in the supplementary material, where we design three types of experiments to verify the basic idea with simulated date.) B. Ionosphere Three-Dimensional Velocity Image Currently, whether 2-D ionosphere models or 3-D ionosphere tomography images, it is assumed that the ionosphere electron density does not change in a certain time span (e.g., 2 h); this is actually unreasonable. On one hand, ignoring the ionosphere change in 2 h will affect the inversion accuracy; on the other hand, the ionosphere velocity image not only can provide more reliable gradient information to other joint-inversion algorithm to improve the accuracy, but it also has important implications in ionosphere electron density forecast. In order to introduce the ionosphere velocity parameters, first, we have to analyze the electron density change characteristics in the regular time span (2 h). We calculate the ionosphere velocity value using IRI 2012 model among 35◦ N–65◦ N, 0◦ E–25◦ E from March 1, 2013 to March 1, 2013 and then subtract the velocity values to statistic the average and standard deviation to investigate the ionosphere changes characteristics in the latitude, longitude, and altitude, respectively. Fig. 4 gives out the comparison. Overall, the velocities surely exist and match each other well before and after 1 h. From the three altitudes in Fig. 4(a)–(c), the larger differences concentrate in 200–400 km; the velocities
in other layers maintained a good compliance. Fig. 4(a) and (b) gives the comparison of ionosphere velocity on 300 and 400 km, respectively; from the figure, we can see that the big differences exist on March 1–3, 2013, with the Kp indexes all beyond 15, but it can ensure that the ionosphere velocities before and after 1 h are on the same order of magnitude in the whole area; thus, it guarantees that using the first-order velocities term in 2 h to estimation is reasonable. In addition, the electron density values are an order of magnitude larger than the velocities term, considering that the first order of velocities does not affect least squares estimation. The preceding discussion shows that it is reasonable to add the ionosphere velocity parameters on the basis of the original ionosphere electron density parameters. However, after considering the grids’ velocity parameters xv , we should take further consideration about the influence on forming observations equation in a certain modeling time span, i.e., Li = STECi (t) =
J
→ (xj + xvj )bj (− r ) ds
s j=1
=
J
xj
j=1
→ bj ( − r ) ds +
J j=1
s
xvj
→ bj ( − r ) ds,
s
(i = 1, 2, . . . , M ) (9) ⎞ ⎛ − → b1(→ r ) ds 0 (t−t0) bi (− r ) ds 0 s . ⎟ ⎜s . . . ⎟ ⎜ ⎟ ⎜ . . . . Li =⎜ ⎟ ⎟ ⎜ . . . . − ⎠ ⎝ → bi (→ r ) ds . (t−t0) bn (− r ) ds . s
⎛
s
⎞
x1 ⎜ −− · ⎟ ⎜ ⎟ ⎜x ⎟ x ⎜ ⎟ n ×⎜ =A = AX. ⎟ x v ⎜ xv1 ⎟ ⎝ . ⎠
(10)
xvn However, the additional velocity parameters increase the condition number of normal equation further and make the estimation more unstable. To overcome this problem, we use the same grids constraints method introduced in Section III-A to add constraint to the grid velocity parameters. First, the initial value of the grid velocity parameters N eV can be obtained by the difference of the grid electron density values, i.e., N eV = N et − N et−1 .
(11)
With the initial velocity parameters, we use the conditions (12) and (13) to add constraint between the grids, i.e., ⎛
⎞ N eV 1 ⎜ N eV 2 ⎟ ⎜ ⎟ Bxv = 0 xv = ⎜ . ⎟ ⎝ ⎠ . N eV n
B=
B1 B2
(12)
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
2377
Fig. 4. Comparison is between 13:00–14:00 and 14:00–15:00 on March 1, 2013, with Sum[Kp] = 35, 13:00–14:00 and 14:00–15:00 on March 2, 2013, with Sum[Kp] = 22+, 13:00–14:00 and 14:00–15:00 on March 3, 2013, with Sum[Kp] = 17+, 08:00–09:00 and 09:00–10:00 on March 4, 2013, with Sum[Kp] = 9−, 21:00–22:00 and 22:00–23:00 on March 5, 2013, with Sum[Kp] = 7−, 17:00–18:00 and 18:00–19:00 on March 6, 2013, with Sum[Kp] = 5+, 21:00–22:00 and 22:00–23:00 on March 7, 2013, with Sum[Kp] = 4+. (a) Comparison of ionosphere velocity on altitude 300 km. (b) Comparison of ionosphere velocity on altitude 400 km. (c) Comparison of ionosphere velocity on altitude 800 km. (d) Comparison of ionosphere velocity on Lat = 35◦ N surface. (e) Comparison of ionosphere velocity on Lat = 55◦ N surface. (f) Comparison of ionosphere velocity on Lat = 65◦ N surface.
where B is the condition matrix, B1 is the condition matrix on the same altitude, and B2 is the condition matrix at the same site. B1 and B2 can be expressed with formula (13), i.e., ⎞ ⎛ I G1 I G2 0 . . 0 . . . . ⎟ ⎜ . ⎟ ⎜ . I Gk 0 . 0⎟ ⎜IG1 B1=⎜ ⎟ . . IGk+1 IGk+2 0 ⎟ ⎜0 ⎠ ⎝ . . . . . . . I Gn 0 . 0 IGk+1 ⎛ I G1 ⎜IG1 ⎜ ⎜0 B2=⎜ ⎜0 ⎝ . 0
. IGk+1 0 . . . . . IG2k+1 . . . . . I G2 . . . 0 I G2 . . . . . . . . 0 .
(13a)
⎞ . . 0 . . 0 ⎟ ⎟ . IGk+1 0 ⎟ ⎟. . . IG2k+2⎟ ⎠ . . . . . . (13b)
C. Anomaly Detection and Repair of Electron Density Estimation Due to the ill-posed problem of the normal equation, there will be a larger bias between the electron density inversion estimation and its real value. In order to prevent the introduction of these outliers as the initial values to the next iteration and then affect the accuracy of the final result, it is necessary to conduct the detection and repair of the abnormal electron density estimation. There are two main aspects for the detection and repair. 1) Anomaly Estimation Detection and Repair on the Same Altitude: To the electron density estimation among different layers, we use a polynomial fitting model to detect and repair the anomaly estimation at the same altitude, i.e., x = a0 + a 1 b + a 2 l + · · · .
(14)
After fitting the electron density estimation, we calculate residuals and modify the electron density values to the obvious abnormalities.
2378
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 5. Flowchart of the new ionosphere tomography algorithm with loose grid constraints and 3-D velocity profile. The left figure gives the data process flow, and the right figure gives the formula process flow.
Fig. 6. Time sequences of the Kp, Dst, and AE indexes from August 14, 2003 to August 15, 2003.
2) Anomaly Estimation Detection and Repair on Different Altitudes at the Same Site: To the electron density estimation among different altitudes at the same site, we use the empirical model Chapman function to do the anomaly detection and repair on different altitudes at the same site, particularly at the edge areas, i.e., N e(h) = N e0 · exp (1 − z − exp(−z)) z = (h − h0)/H
(15)
where N e0 is the electron density peak value, whereas h0 is the corresponding peak height. H is the scale height. After the above two steps repair work, it can basically detect and repair the anomaly estimation and then provide the more reliable initial values for the next iteration and improve the inversion accuracy. The parameter estimation process is illustrated in Fig. 5. IV. N UMERICAL E XPERIMENTS The geomagnetic condition is relatively quiet on August 15, 2003. Fig. 6 gives the time sequence of the Kp, Dst, and Auroral Electrojet (AE) indexes from August 14, 2003 to August 15, 2003. The Kp index is below 4, and the Dst index is lower than −20 nT. The AE index is almost below 300 on August 15. In this paper, we use the European GNSS data of the 38 IGS stations and reconstruct the ionosphere electron density and
Fig. 7.
Schematic of the used IGS station distribution.
velocity 3-D images from 0◦ E–20◦ E, 35◦ N–55◦ N. Fig. 7 gives the station distribution of the used station, whereas Fig. 8 gives the variations images of the electron density and velocity on longitude 12◦ surfaces. Fig. 8 basically reflects the ionosphere changes of a day under quiet geomagnetic condition. From the images, it can be seen that the peak height of electron density experiences a
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
2379
Fig. 8. Electron density and velocity images on 12◦ E surface on August 15, 2003 (unit: Ne/m3 ). The left images are the electron density, and the right images are the electron density velocity images.
2380
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 8. (Continued) Electron density and velocity images on 12◦ E surface on August 15, 2003 (unit: Ne/m3 ). The left images are the electron density, and the right images are the electron density velocity images.
decline and then a rise process. The process reflects the diurnal variation of ionosphere electron density vertical structure. As can be seen from the velocity images, the main electron density activities focused on 200- to 500-km range. Here, the ionosphere change of a day can be discussed in four stages.
TABLE I R ESULTS ’ C OMPARISON A MONG D IFFERENT L ATITUDE S URFACES B ETWEEN D IFFERENT M ETHODS (U NIT: 1010 Ne/m3 ) M ODELING T IME: 13:00–15:00
• The primary change phase: The quiet stage (00:00UT– 06:00UT) In Fig. 8(a) specifically, the peak height of electron density is at 350-km surface height at UT02:00, and the electron density peak value is smaller than 4.0 × 1011 el/m3 . The velocity image in Fig. 8(b) shows that the electron density starts to change while the velocity of 200–300 km increases and the velocity of 300- to 500-km height decreases, but the magnitude is smaller than 1.0 × 1011 el/m3 . Thus, the ionosphere is relatively calm at this stage. • The second change phase: The climbing stage (06:00UT– 12:00UT) Driven by the primary phase’s influence, the electron density peak height decreased to 250-km surface at UT06:00 in Fig. 8(c), the electron density magnitude increases at this time, and the velocity image in Fig. 8(d) has similar structural characteristics with the electron density image. The velocity variation rate significantly increases comparing with UT02:00. Overall, the electron density experiences the increase process;
in particular, the velocity at 200 km reaches 1.0 × 1011 el/m3 magnitude between 35N and 40N. The electron density continues to rise at UT10:00 in Fig. 8(e), the peak height experiences the rebound phenomena. In Fig. 8(f), the velocity at 200 km is converted to be negative, and the positive velocity at 300 km tends to be strong.
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
Fig. 9.
2381
Electron density change sequence among different layers around UT18:00 of ionosonde, IRI, and UpIRI data.
• The third phase: The maximum stage (12:00UT–15:00UT) At UT14:00 in Fig. 8(g), the electron density peak reaches the maximum magnitude of the day at 1.0 × 1012 el/m3 , and the peak height rises back to 300 km. Above 300 km height, the electron density change is smooth. In addition, the electron density reaches maximum, but it needs to be noted that the electron density decreases in greater strength below 200 km from 35N to 55N now. • The fourth phase: The declining stage and quiet stage (16:00UT–24:00UT) Comparing with UT14:00, the electron density magnitude decreases, but the peak height keeps at 300 km at UT18:00 in Fig. 8(i). There is a negative velocity at peak height. From the velocity image, another point worth noting is that there is a positive velocity at 4.0 × 1010 el/m3 magnitude at 350 km from 45N to 55N; driven by this, at UT22:00, in the case of the electron density reduction in other inversion area, the electron density at 350 km from 45N to 55N increases comparing with the UT18:00 image. This distinctive increase phenomenon proves the feasibility of forecasting the 3-D ionosphere change with 3-D velocity images. We will verify this phenomenon with ionosonde data next. At UT22:00, the electron density peak goes back to the 350-km height, from Fig. 8(l), at 250–350 km, the electron density decreases in maximum strength of the day, but the changes of other ionosphere layers are relatively smooth. On the whole, taking 4:00 hour’s time resolution into account, the electron density variation is in good agreement with the velocity images, and it shows that this study presents a
way to forecast the 3-D electron density variation with the 3-D velocity images. Furthermore, we use the conventional algorithm MART to reconstruct the ionosphere electron density images at the same area and the same time and compare the MART’s results with the new algorithm’s results among different latitude surfaces. The comparisons are shown in Table I. In the table, it is shown that with the latitudes reduced, the differences of maximum value, averages, and standard deviations are gradually increasing; specifically, the maximum value is increased by one order of magnitude. However, considering Fig. 8(g), the electron density peak value reaches 1012 el/m3 at 14:00UT, while the differences are much less than the magnitude. That means that the inversion results of the two algorithms are basically consistent. In addition, the ionosphere ionosonde data are downloaded and used to prove the new algorithm’s results. As aforementioned, the new algorithm detects that, on 12◦ E surface at UT18:00, between 45◦ N and 55◦ N, there is an electron density increase phenomenon at 300–400 km, which is inconsistent with the overall change. This paper selects ionosonde data at European regions to verify the reliability of this phenomenon. In this area, there are three ionosonde observatory stations, but only JR055 JULIUSRUH (54.60◦ N, 13.40◦ E) conducted ionosphere exploration that day. In Fig. 8(j), the electron density increases at 300–400 km and decreases at 100–200 km, respectively. JR055 conducts 11 measurements from UT17:28 to UT20:13. We intercept 100- to 200-km and 300- to 400-km detection results in Fig. 9. In the figure, the IRI model’s results and the new algorithm’s results (the UpIRI) are also printed out as the reference.
2382
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 10. Electron density profiles of different models at 54.60◦ N, 13.40◦ E. The top panel is the electron density profile at UT14:00, whereas the bottom panel is at UT20:00.
As shown in the figure, the detection results of ionosonde have exactly the same variation process comparing with the new algorithm’s results. Particularly on layers 300–400 km, the IRI model fails to reflect the real electron density increasing phenomena, whereas the new algorithm properly reflects the abnormal disturbance with high electron density detection accuracy. The velocity images of the new algorithm on UT18:00 are consistent with the ionosonde’s detection results. Fig. 10 shows the electron density profiles of different models at the JR055 station site. In the top panel in Fig. 10, the electron density peak height experiences a decline, comparing the ionosonde data with the IRI model (from 350 to 300 km). The new algorithm result shows the peak height decline exactly. The bottom panel gives the comparison at UT20:00. The electron density value of ionosonde data is greatly increased comparing with the IRI model on layers 200–500 km, whereas the new algorithm gives more precise results. V. C ONCLUSION The real-time detection and analysis of the Earth’s space environment abnormal outbreak activity caused by some natural disasters are likely to be a potential and ancillary means to forecast major natural disasters. The ionosphere tomography technology based on GNSS can obtain the 3-D ionosphere morphology. The inversion algorithm is the critical challenge, which restricts the precision and accuracy of the ionosphere tomography technology. In this paper, by adding the loose constraints among the grids,
forming a stable mesh structure and then composing the virtual normal equation, combining the virtual normal equation with the observed normal equation, the electron density parameters are estimated with best fit to the real observations based on the priori information. In addition, different from previous algorithms, the new algorithm introduces the electron density velocity parameters, which is more reasonable with actual ionosphere variation. The 3-D velocity profiles are useful and provide multisource information for studying the spatial ionosphere variation, particularly the ionosphere variation under ionosphere disturbances. In addition, the 3-D ionosphere changes velocity can provide reliable ionosphere gradient information in each layer for some joint-inversion algorithms and then improve the accuracy of ionosphere model. We have used the new algorithm to reconstruct the ionosphere electron density and velocity images on August 15, 2003, around Europe. By comparing the results with the conventional algorithm MART’s results and ionosphere ionosonde data, the feasibility and reliability of the new algorithm are verified. The related research provides a new way to real-time detection and prediction of ionosphere change. ACKNOWLEDGMENT The authors would like to thank IGS for providing GNSS data and global ionospheric radio observatory (GIRO) for providing the ionosonde data. The authors would also like to thank the anonymous reviewers for their constructive comments in improving this paper.
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
R EFERENCES [1] M. M. Alizadeh, H. Schuh, S. Todorova, and M. Schmidt, “Global ionosphere maps of VTEC from GNSS, satellite altimetry, formosat3/COSMIC data,” J. Geodesy, vol. 85, no. 12, pp. 975–987, Dec. 2011. [2] C. K. Avinash and S. Malcolm, Principles of Computerized Tomographic Imaging. Piscataway, NJ, USA: IEEE Press, 1988, pp. 275–296. [3] K. Davies and D. M. Baker, “Ionospheric effect observed around the time of the Alaskan earthquake of March 28, 1964,” J. Geophys. Res., vol. 70, no. 9, pp. 2251–2253, May 1965. [4] R. S. Leonard and R. A. Barnes, “Observation of ionospheric disturbances following the Alaska earthquake,” J. Geophys. Res., vol. 70, no. 5, pp. 1250–1253, Mar. 1965. [5] D. Dettmering, M. Schmidt, and R. Heinkelmann, “Combination of different space-geodetic observations for regional ionosphere modeling,” J. Geodesy, vol. 85, no. 12, pp. 989–998, Dec. 2011. [6] R. Garcia and F. Crespon, “Radio tomography of the ionosphere: Analysis of an underdetermined, ill-posed inverse problem, regional application,” Radio Sci., vol. 43, no. 2, pp. RS2014-1–RS2014-13, Apr. 2008. [7] T. Hobiger, T. Kondo, and Y. Koyama, “Constrained simultaneous algebraic reconstruction technique (C-SART)—A new and simple algorithm applied to ionospheric tomography,” Earth Planets Space, vol. 60, no. 7, pp. 727–735, Jul. 2008. [8] Z. Z. Liu, “Ionosphere tomographic modeling and applications using Global Positioning System (GPS) measurements,” Univ. Calgary, Calgary, AB, USA, 2004. [9] J. Y. Liu et al., “Pre-earthquake ionospheric anomalies registered by continuous GPS TEC measurements,” Ann. Geophys.-Germany, vol. 22, no. 5, pp. 1585–1593, 2004. [10] E. P. Macalalad, L. C. Tsai, J. Wu, and C. H. Liu, “Application of the Taiwan Ionospheric Model to single-frequency ionospheric delay corrections for GPS positioning,” GPS Solutions, vol. 17, no. 3, pp. 337–346, Jul. 2013. [11] S. Kon, M. Nishihashi, and K. Hattori, “Ionospheric anomalies possibly associated with M ≥ 6.0 earthquakes in the Japan area during 1998–2010: Case studies and statistical study,” J. Asian Earth Sci., vol. 41, no. 4/5, pp. 410–420, Jun. 2011. [12] A. Krankowski, I. Zakharenkova, A. Krypiak-Gregorczyk, I. I. Shagimuratov, and P. Wielgosz, “Ionospheric electron density observed by FORMOSAT-3/COSMIC over the European region and validated by ionosonde data,” J. Geodesy, vol. 85, no. 12, pp. 949–964, Dec. 2011. [13] P. F. Fougere, “Ionospheric radio tomography using maximum entropy,” Radio Sci., vol. 30, no. 2, pp. 429–444, Mar./Apr. 1995. [14] S. E. Pryse and L. Kersley, “A preliminary experimental test of ionospheric tomography,” J. Atmos. Terrestrial Phys., vol. 54, no. 7/8, pp. 1007–1012, Jul./Aug. 1992. [15] S. A. Pulinets and K. Boyarchuk, Ionospheric Precursors of Earthquakes. Berlin, Germany: Springer-Verlag, 2004. [16] M. Schmidt, D. Bilitza, C. K. Shum, and C. Zeilhofer, “Regional 4D modeling of the ionospheric electron density,” Adv. Space Res., vol. 42, pp. 782–790, 2008. [17] S. Saroso, J. Y. Liu, K. Hattori, and C. H. Chen, “Ionospheric GPS TEC anomalies and M > 5.9 earthquakes in Indonesia during 1993–2002,” Terrestrial Atmos. Ocean. Sci., vol. 19, no. 5, pp. 481–488, Oct. 2008. [18] S. Todorova, H. Schuh, and T. Hobiger, “Using the global navigation satellite systems and satellite altimetry for combined global ionosphere maps,” Adv. Space Res., vol. 42, no. 4, pp. 727–736, Aug. 2007. [19] L. C. Tsai, W. H. Tsai, W. S. Schreiner, F. T. Berkey, and J. Y. Liu, “Comparisons of GPS/MET retrieved ionospheric electron density and ground based ionosonde data,” Earth Planets Space, vol. 53, pp. 193–205, 2001.
2383
[20] D. B. Wen and Y. B. Yuan, “A hybrid reconstruction algorithm for 3-D ionospheric tomography,” IEEE Trans.Geosci. Remote Sens., vol. 46, no. 6, pp. 1733–1739, Jun. 2008. [21] D. B. Wen and Y. Wang, “A new two-step algorithm for ionospheric tomography solution,” GPS Solutions, vol. 16, no. 1, pp. 89–94, Jan. 2012. [22] Y. B. Yao et al., “Analysis of pre-earthquake ionospheric anomalies before the global M = 7.0+ earthquakes in 2010,” Nat. Hazards Earth Syst. Sci., vol. 12, no. 3, pp. 575–585, 2012. [23] X. N. Yue, W. S. Schreiner, and Y. H. Kuo, “Evaluating the effect of the global ionospheric map on aiding retrieval of radio occultation electron density profiles,” GPS Solutions, vol. 17, no. 3, pp. 327–335, Jul. 2013.
Yibin Yao received the B.Sc., Master’s, and Ph.D. degrees (with distinction) in geodesy and surveying engineering from Wuhan University, Wuhan, China, in 1997, 2000, and 2004, respectively. He is currently a Professor with Wuhan University. His main research interests include global navigation satellite system ionospheric/atmospheric/ meteorological studies, theory and method of surveying data processing, and GPS/MET and highprecision GPS data processing. He has authored over 50 publications since 2001 and six books in those areas. Prof. Yao has been a recipient of a Luojia Professorship at Wuhan University, supported by the Program for New Century Excellent Talents in University.
Jian Kong received the B.Sc. degree from Shandong University of Science and Technology, Qingdao, China, in 2009. He is currently working toward the Ph.D. degree in geodesy and surveying engineering in the School of Geodesy and Geomatics, Wuhan University, Wuhan, China. His main research interests include GNSS data processing, global ionospheric modeling using multisource geodesy observations, 3-D ionospheric tomography based on GNSS observations, and ionospheric anomalies analysis under abnormal conditions.
Jun Tang received the B.Sc. degree in surveying and mapping engineering and the M.Sc. degree in geodesy and surveying engineering from East China Institute of Technology, Fuzhou, China, in 2006 and 2009, respectively. He is currently working toward the Ph.D. degree in geodesy and surveying engineering in the School of Geodesy and Geomatics, Wuhan University, Wuhan, China. His research interests include ionospheric modeling using diverse data, multidimensional ionospheric tomography based on global navigation satellite system (GNSS) observations, and GNSS applications.
2373
A New Ionosphere Tomography Algorithm With Two-Grid Virtual Observations Constraints and Three-Dimensional Velocity Profile Yibin Yao, Jian Kong, and Jun Tang
Abstract—Ionosphere tomography is a typical ill-posed problem, and using the ionosphere priori information as the constraints to improve the state of normal equation is an effective approach to solve this problem. In this paper, we impose priori constraints by increasing the virtual observations in n-dimensional space. Then, after the inversion region to be gridded, we can form a stable structure between the grids with loose constraints, which greatly improve the state of normal equation. Based on that, to obtain the real-time velocity information of ionosphere electron density, we introduce the grid electron density velocity parameters, which can be estimated with electron density parameters simultaneously. The authors use the new algorithm and the global navigation satellite system data in Europe to inverse the 12 ionosphere electron density and velocity images on August 15, 2003, which reflects the ionosphere changes of electron density and velocity in whole day; in addition, we compare the results with the traditional algorithm multiplicative algorithm reconstruction technique’s results. Furthermore, we compare the results with the changes time series of plasma frequency observed by ionosphere ionosonde among different layers, and many types of analysis and comparison verify the effectiveness and reliability of the new algorithm. The related research provides a new way for the real-time detection and prediction of ionosphere changes. Index Terms—Computerized ionosphere tomography, grid constraints, three-dimensional ionosphere velocity image, virtual observations.
I. I NTRODUCTION
A
S THE important part of the Earth space environment, ionosphere about 60–1000 km above the ground is the critical protective layer for human survival and development. Using different kinds of modern technologies to detect and analyze environment changes of the Sun–Earth space, particularly the Earth’s space, and then to investigate their basic structure and variation helps us not only to improve accuracy of the velocimetry, location, timing, precision communications,
Manuscript received April 17, 2014; revised July 27, 2014; accepted September 8, 2014. This work was supported in part by the National Natural Science Fund of China under Grant 41274022 and Grant 41174012 and in part by the 863 Program under Grant 2013AA122502. This paper has supplementary downloadable material available at http://ieeexplore.ieee.org, provided by the authors. This includes details about the basic principle of the new algorithm and verification of the basic idea with simulated data. This material is 214 kB in size. The authors are with the School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2014.2359762
and navigation but also to study and find out the relationship between the upper atmospheric layers. Particularly, it has important scientific significance in the study of the mechanism of global ionosphere disturbances and irregular changes. Ionosphere tomography is a typical ill-posed problem. Such problems as fewer observations in the edge of the inversion region and the observations from global navigation satellite system (GNSS) constellation are nearly vertical lead to the condition number of the normal equation in inversion algorithm being too large, and then, the estimation solution is unstable with the presence of observation noise. The commonly used reconstruction algorithms are algorithm reconstruction technique (ART), multiplicative ART (MART), and simultaneous iteration reconstruction technique (SIRT) [2], [13], [14]. The most intuitive idea to solve the ill-posed problem is increasing the observations, such as, currently, the ionosphere inversion algorithms based on multisource data [1], [5], [12], [16], [18], [19]. On the other hand, introducing the priori information of ionosphere electron density to improve the state of normal equation also can overcome the ill-posed problem. Hobiger et al. [7] used the priori constraint with SIRT to inverse ionosphere electron density and then compared the results with SIRT. After imposing the priori constraints, the iterative convergence speed and accuracy are improved. In contrast to the lack of constraint in the edge area, Liu [8] proposed an improved SIRT algorithm with priori constraints and verified its superiority. In recent years, Garcia and Crespon, Hobiger et al., Wen and Yuan, and Wen and Wang have proposed several corresponding solutions to ill-posed problems in ionosphere tomography due to insufficient projection data [6], [7], [20], [21]. However, it needs to be noted that the ill-posed problem of normal equation is still the critical problem in ionosphere tomography algorithm. Ionosphere anomaly detection, particularly the related ionosphere anomalies and ionosphere scintillation before the natural disasters (tsunamis, earthquakes), is a more popular research direction [15]. However, most of the research are confined in studying the statistical analysis of the spatial location and time position of the anomalies, which is based on the 2-D ionosphere images [3], [4], [8], [11], [17], [22]; there are fewer research about the physical mechanisms and theoretical models of the anomalies based on 3-D ionosphere images. The establishment of physical model needs the reliable and adequate ionosphere change information. The existing information for model building is insufficient to provide further information,
0196-2892 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
2374
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 1. Schematic of the computation domain of tomography.
and ionosphere electron density information in real time (ionosphere changes velocity) has an important role in the physical model establishment. With the 3-D ionosphere electron density images, the 3-D ionosphere changes velocity images can provide reliable ionosphere gradient information in each layers for some joint inversion algorithms and then improve the accuracy of ionosphere model and positioning [13], [23]. On the other hand, the 3-D ionosphere changes velocity information is conducive to ionosphere prediction and forecasting in three dimensions, which is of great significance to scientific research and engineering applications based on ionosphere changes. II. BASIC P RINCIPLES OF I ONOSPHERE E LECTRON D ENSITY R ECONSTRUCTION BASED ON THE P IXEL L EVEL The electron density reconstruction is a typical inverse problem, which is based on the observed slant total electron content (STEC) to reconstruct the 3-D distribution in formation of ionosphere electron density. As shown in Fig. 1, after the inversion region to be gridded in three dimensions, GPS satellite signal propagation paths can be taken approximately as a straight line. → Select an appropriate set of basic functions b(− r ) to model the ionosphere electron density in the inversion region, i.e., → N e(− r , t) ≈
J
→ xj (t) · bj (− r)
(1)
j=1
where J represents the number of selected basis functions, xj (j = 1, 2, . . . , J) represents the coefficients of the basis → function, N e(− r , t) is the ionosphere electron density on the → signal propagation paths, t is the observation time, and − r is the ray path vector. Then, the STEC of each ray path can be expressed as follows: STECi (t) =
J
→ xj (t)bj (− r ) ds
s j=1
=
J j=1
xj (t) s
→ bj ( − r ) ds,
i = 1, 2, . . . , M
(2)
where M represents STEC number (the GPS ray path num the → r ) ds, considering the GPS measurement ber). Let Aij = l bj (− noise, then the STEC on GPS signal propagation paths can be expressed as follows: STECi =
J
Aij xj + εi .
(3)
j=1
The formula (3) can be expressed by the vector equation as in formula (4), i.e., L = STEC = Ax · x + ε
(4)
where L is the observation vector. STEC is the column vector of STEC. Ax is the design matrix. x is the column vector composed of the coefficients of basic function. In the ionosphere tomography, x is the ionosphere electron density. ε is the column vector of the observation noise. III. N EW I ONOSPHERE T OMOGRAPHY A LGORITHM W ITH T WO -G RID C ONSTRAINTS AND T HREE -D IMENSIONAL V ELOCITY P ROFILE A. Loose Constraints Among Grids In the process of GNSS-based electron density inversion in pixel level, on one hand, due to the lack of observational data, there are no rays that pass through some grids, and on the other hand, the GNSS observations rays are almost vertical; there is limited horizontal information carried by the observation rays. Affected by the two reasons, the condition number of the normal equation matrix in parameter estimation is very large; thus, the inversion of matrix is unstable and inaccurate. There are two main ways to solve this problem. First, with the maturity and development of various ionosphere detection satellites (DORIS, occultation, satellite altimetry), the observational data types that can be used are gradually increasing. The joint inversion combining multisource data is an effective way to solve the preceding problems. In addition, although the ionosphere is affected by the Sun, geomagnetic activity, and other factors, as the important part of Sun–Earth environment, the ionosphere changes reflect the electron content change, which is the electron moving and transmission process. That
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
2375
Fig. 3. Schematic of plane constraint and altitude constraint.
Fig. 2.
Cooperation of the two different ways to add virtual observations.
means that the electron density values among different grids are correlated with each other in time and space scales. Thus, based on this, adding the constraints among the grids and then improving the state of normal equation is another effective solution. There are many ways to add constraints, such as increasing parameters constraint equations. However, the more simple and practical method is to add constraints by increasing the virtual observation to form the virtual normal equation. According to the external auxiliary information [such as global ionosphere map (GIM) and international reference ionosphere (IRI)], simulating virtual station location, virtual satellite position, and virtual observations and then combining the virtual and real observations to consist the normal equations simultaneously improves the state of the equation. Increasing the virtual observations can solve the problem of insufficient observed data in some extent, but it still cannot solve the problem that the observation data are nearly vertical, because the virtual observation values are still lacking the constraint in the horizontal direction. Thus, after adding the virtual observation, the condition number of normal equation is still large. Based on the existing algorithm, in this paper, we propose a new way to add the virtual observations, in order to increase the horizontal constraint in the normal equation and improve its state. The traditional method [see Fig. 2(a)] to add virtual observation is to simulate station location, satellite positions, and STEC; and this method is implemented in 3-D space. If we do it again in n-dimensional space, a ray from the satellite can pass through any two arbitrary grids to reach stations in
n-dimensional space [see Fig. 2(b)]. This new way to add virtual observations cannot only ensure that each grid can be passed through by virtual rays but also increase the horizontal constraint among the grids by designing the virtual rays. Then, combine the virtual observations in n dimension and the real observations to form the normal equation for parameter estimation. Compared with the traditional method, the state improvement to the normal equation of the new approach is very clear. However, if we simulate a virtual ray between every two grids, the constraints among the grids will be too tight; in case the real observations are insufficient, the grid system with tight constraints cannot be corrected rationally to keep with the real observations. Considering that the ionosphere inversion algorithms are mostly based on altitude and horizontal direction, in this paper, we add virtual observations on these two levels based on this characteristic, which is shown in Fig. 3. • horizontal constraint C: add the plane constraints in horizontal direction; • altitudinal constraint H: add the proportion constraints in altitude direction where H and C are determined by the initial value of the grid electron density. To add the H and C constraints to the estimation, first, select the grids Gi and Gj on the same altitude or at the same site orderly, and then, obtain the two grids’ electron densities N eGi and N eGj according to the IRI model. Based on this, simulate two intercepts IGi and IGj crossing the two grids artificially; here, we need to refer to the real grid size and the size range of the actual intercepts. Then, a virtual observation ¯l can be obtained with the formula (5), i.e., ¯l = IG • N eG + IGj • N eGj . i i
(5)
After forming the virtual observations among all the two grids on the same surfaces and the same sites with formula (5), ¯ i.e., we can get the virtual observations matrix L, ⎛¯ ⎞⎛ ⎞ I G1 • N e G1 + I G2 • N e G2 l1 . ⎟⎜ . ⎟ ¯=⎜ L (6) ⎝ ⎠⎝ ⎠. . . ¯ln IGn−1 • N eGn−1 + IGn • N eGn
2376
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Then, form the virtual observations to the normal equation, i.e., ⎞⎛ ⎞ ⎛ x1 0 0 IG1 IG2 . ⎟ . . . . ⎟ ⎜ ⎜ ¯=⎝ L ⎠⎝ ⎠ . . . . . xn 0 0 IGn−1 IGn ⎛ ¯ A11 ⎜ . =⎝ . 0
A¯12 . . 0
0 . . A¯nn−1
⎞⎛ ⎞ x1 0 . ⎟⎜ . ⎟ ⎠ ⎝ ⎠ = A¯x x. (7) . . A¯nn xn
The simultaneous equations can be obtained with formula (4) and formula (7), i.e., Ax L PL 0 (8) ¯ = A¯x x P = L 0 PL¯ where P is the weight matrix, PL is the weight matrix for the real observations, and PL¯ is the weight matrix for the virtual observations. In this paper, we set PL /PL¯ = 100/1. After adding these two constraints, the inversion area is formed as a stable structure with loose constraints among grids. Then, if there is a real ray passing through the area, with greater weight on the real observation, then the grid the rays pass through will selfcorrect to match real observation. Due to the plane and altitude constraints among the grids, the corrected grids will trigger the other grids in the inversion area to self-correction order to comply optimally with the real observation data; eventually, the whole inversion region achieves optimal correction based on actual observations. (More details about the basic principle of the new algorithm can be found in the supplementary material, where we design three types of experiments to verify the basic idea with simulated date.) B. Ionosphere Three-Dimensional Velocity Image Currently, whether 2-D ionosphere models or 3-D ionosphere tomography images, it is assumed that the ionosphere electron density does not change in a certain time span (e.g., 2 h); this is actually unreasonable. On one hand, ignoring the ionosphere change in 2 h will affect the inversion accuracy; on the other hand, the ionosphere velocity image not only can provide more reliable gradient information to other joint-inversion algorithm to improve the accuracy, but it also has important implications in ionosphere electron density forecast. In order to introduce the ionosphere velocity parameters, first, we have to analyze the electron density change characteristics in the regular time span (2 h). We calculate the ionosphere velocity value using IRI 2012 model among 35◦ N–65◦ N, 0◦ E–25◦ E from March 1, 2013 to March 1, 2013 and then subtract the velocity values to statistic the average and standard deviation to investigate the ionosphere changes characteristics in the latitude, longitude, and altitude, respectively. Fig. 4 gives out the comparison. Overall, the velocities surely exist and match each other well before and after 1 h. From the three altitudes in Fig. 4(a)–(c), the larger differences concentrate in 200–400 km; the velocities
in other layers maintained a good compliance. Fig. 4(a) and (b) gives the comparison of ionosphere velocity on 300 and 400 km, respectively; from the figure, we can see that the big differences exist on March 1–3, 2013, with the Kp indexes all beyond 15, but it can ensure that the ionosphere velocities before and after 1 h are on the same order of magnitude in the whole area; thus, it guarantees that using the first-order velocities term in 2 h to estimation is reasonable. In addition, the electron density values are an order of magnitude larger than the velocities term, considering that the first order of velocities does not affect least squares estimation. The preceding discussion shows that it is reasonable to add the ionosphere velocity parameters on the basis of the original ionosphere electron density parameters. However, after considering the grids’ velocity parameters xv , we should take further consideration about the influence on forming observations equation in a certain modeling time span, i.e., Li = STECi (t) =
J
→ (xj + xvj )bj (− r ) ds
s j=1
=
J
xj
j=1
→ bj ( − r ) ds +
J j=1
s
xvj
→ bj ( − r ) ds,
s
(i = 1, 2, . . . , M ) (9) ⎞ ⎛ − → b1(→ r ) ds 0 (t−t0) bi (− r ) ds 0 s . ⎟ ⎜s . . . ⎟ ⎜ ⎟ ⎜ . . . . Li =⎜ ⎟ ⎟ ⎜ . . . . − ⎠ ⎝ → bi (→ r ) ds . (t−t0) bn (− r ) ds . s
⎛
s
⎞
x1 ⎜ −− · ⎟ ⎜ ⎟ ⎜x ⎟ x ⎜ ⎟ n ×⎜ =A = AX. ⎟ x v ⎜ xv1 ⎟ ⎝ . ⎠
(10)
xvn However, the additional velocity parameters increase the condition number of normal equation further and make the estimation more unstable. To overcome this problem, we use the same grids constraints method introduced in Section III-A to add constraint to the grid velocity parameters. First, the initial value of the grid velocity parameters N eV can be obtained by the difference of the grid electron density values, i.e., N eV = N et − N et−1 .
(11)
With the initial velocity parameters, we use the conditions (12) and (13) to add constraint between the grids, i.e., ⎛
⎞ N eV 1 ⎜ N eV 2 ⎟ ⎜ ⎟ Bxv = 0 xv = ⎜ . ⎟ ⎝ ⎠ . N eV n
B=
B1 B2
(12)
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
2377
Fig. 4. Comparison is between 13:00–14:00 and 14:00–15:00 on March 1, 2013, with Sum[Kp] = 35, 13:00–14:00 and 14:00–15:00 on March 2, 2013, with Sum[Kp] = 22+, 13:00–14:00 and 14:00–15:00 on March 3, 2013, with Sum[Kp] = 17+, 08:00–09:00 and 09:00–10:00 on March 4, 2013, with Sum[Kp] = 9−, 21:00–22:00 and 22:00–23:00 on March 5, 2013, with Sum[Kp] = 7−, 17:00–18:00 and 18:00–19:00 on March 6, 2013, with Sum[Kp] = 5+, 21:00–22:00 and 22:00–23:00 on March 7, 2013, with Sum[Kp] = 4+. (a) Comparison of ionosphere velocity on altitude 300 km. (b) Comparison of ionosphere velocity on altitude 400 km. (c) Comparison of ionosphere velocity on altitude 800 km. (d) Comparison of ionosphere velocity on Lat = 35◦ N surface. (e) Comparison of ionosphere velocity on Lat = 55◦ N surface. (f) Comparison of ionosphere velocity on Lat = 65◦ N surface.
where B is the condition matrix, B1 is the condition matrix on the same altitude, and B2 is the condition matrix at the same site. B1 and B2 can be expressed with formula (13), i.e., ⎞ ⎛ I G1 I G2 0 . . 0 . . . . ⎟ ⎜ . ⎟ ⎜ . I Gk 0 . 0⎟ ⎜IG1 B1=⎜ ⎟ . . IGk+1 IGk+2 0 ⎟ ⎜0 ⎠ ⎝ . . . . . . . I Gn 0 . 0 IGk+1 ⎛ I G1 ⎜IG1 ⎜ ⎜0 B2=⎜ ⎜0 ⎝ . 0
. IGk+1 0 . . . . . IG2k+1 . . . . . I G2 . . . 0 I G2 . . . . . . . . 0 .
(13a)
⎞ . . 0 . . 0 ⎟ ⎟ . IGk+1 0 ⎟ ⎟. . . IG2k+2⎟ ⎠ . . . . . . (13b)
C. Anomaly Detection and Repair of Electron Density Estimation Due to the ill-posed problem of the normal equation, there will be a larger bias between the electron density inversion estimation and its real value. In order to prevent the introduction of these outliers as the initial values to the next iteration and then affect the accuracy of the final result, it is necessary to conduct the detection and repair of the abnormal electron density estimation. There are two main aspects for the detection and repair. 1) Anomaly Estimation Detection and Repair on the Same Altitude: To the electron density estimation among different layers, we use a polynomial fitting model to detect and repair the anomaly estimation at the same altitude, i.e., x = a0 + a 1 b + a 2 l + · · · .
(14)
After fitting the electron density estimation, we calculate residuals and modify the electron density values to the obvious abnormalities.
2378
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 5. Flowchart of the new ionosphere tomography algorithm with loose grid constraints and 3-D velocity profile. The left figure gives the data process flow, and the right figure gives the formula process flow.
Fig. 6. Time sequences of the Kp, Dst, and AE indexes from August 14, 2003 to August 15, 2003.
2) Anomaly Estimation Detection and Repair on Different Altitudes at the Same Site: To the electron density estimation among different altitudes at the same site, we use the empirical model Chapman function to do the anomaly detection and repair on different altitudes at the same site, particularly at the edge areas, i.e., N e(h) = N e0 · exp (1 − z − exp(−z)) z = (h − h0)/H
(15)
where N e0 is the electron density peak value, whereas h0 is the corresponding peak height. H is the scale height. After the above two steps repair work, it can basically detect and repair the anomaly estimation and then provide the more reliable initial values for the next iteration and improve the inversion accuracy. The parameter estimation process is illustrated in Fig. 5. IV. N UMERICAL E XPERIMENTS The geomagnetic condition is relatively quiet on August 15, 2003. Fig. 6 gives the time sequence of the Kp, Dst, and Auroral Electrojet (AE) indexes from August 14, 2003 to August 15, 2003. The Kp index is below 4, and the Dst index is lower than −20 nT. The AE index is almost below 300 on August 15. In this paper, we use the European GNSS data of the 38 IGS stations and reconstruct the ionosphere electron density and
Fig. 7.
Schematic of the used IGS station distribution.
velocity 3-D images from 0◦ E–20◦ E, 35◦ N–55◦ N. Fig. 7 gives the station distribution of the used station, whereas Fig. 8 gives the variations images of the electron density and velocity on longitude 12◦ surfaces. Fig. 8 basically reflects the ionosphere changes of a day under quiet geomagnetic condition. From the images, it can be seen that the peak height of electron density experiences a
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
2379
Fig. 8. Electron density and velocity images on 12◦ E surface on August 15, 2003 (unit: Ne/m3 ). The left images are the electron density, and the right images are the electron density velocity images.
2380
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 8. (Continued) Electron density and velocity images on 12◦ E surface on August 15, 2003 (unit: Ne/m3 ). The left images are the electron density, and the right images are the electron density velocity images.
decline and then a rise process. The process reflects the diurnal variation of ionosphere electron density vertical structure. As can be seen from the velocity images, the main electron density activities focused on 200- to 500-km range. Here, the ionosphere change of a day can be discussed in four stages.
TABLE I R ESULTS ’ C OMPARISON A MONG D IFFERENT L ATITUDE S URFACES B ETWEEN D IFFERENT M ETHODS (U NIT: 1010 Ne/m3 ) M ODELING T IME: 13:00–15:00
• The primary change phase: The quiet stage (00:00UT– 06:00UT) In Fig. 8(a) specifically, the peak height of electron density is at 350-km surface height at UT02:00, and the electron density peak value is smaller than 4.0 × 1011 el/m3 . The velocity image in Fig. 8(b) shows that the electron density starts to change while the velocity of 200–300 km increases and the velocity of 300- to 500-km height decreases, but the magnitude is smaller than 1.0 × 1011 el/m3 . Thus, the ionosphere is relatively calm at this stage. • The second change phase: The climbing stage (06:00UT– 12:00UT) Driven by the primary phase’s influence, the electron density peak height decreased to 250-km surface at UT06:00 in Fig. 8(c), the electron density magnitude increases at this time, and the velocity image in Fig. 8(d) has similar structural characteristics with the electron density image. The velocity variation rate significantly increases comparing with UT02:00. Overall, the electron density experiences the increase process;
in particular, the velocity at 200 km reaches 1.0 × 1011 el/m3 magnitude between 35N and 40N. The electron density continues to rise at UT10:00 in Fig. 8(e), the peak height experiences the rebound phenomena. In Fig. 8(f), the velocity at 200 km is converted to be negative, and the positive velocity at 300 km tends to be strong.
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
Fig. 9.
2381
Electron density change sequence among different layers around UT18:00 of ionosonde, IRI, and UpIRI data.
• The third phase: The maximum stage (12:00UT–15:00UT) At UT14:00 in Fig. 8(g), the electron density peak reaches the maximum magnitude of the day at 1.0 × 1012 el/m3 , and the peak height rises back to 300 km. Above 300 km height, the electron density change is smooth. In addition, the electron density reaches maximum, but it needs to be noted that the electron density decreases in greater strength below 200 km from 35N to 55N now. • The fourth phase: The declining stage and quiet stage (16:00UT–24:00UT) Comparing with UT14:00, the electron density magnitude decreases, but the peak height keeps at 300 km at UT18:00 in Fig. 8(i). There is a negative velocity at peak height. From the velocity image, another point worth noting is that there is a positive velocity at 4.0 × 1010 el/m3 magnitude at 350 km from 45N to 55N; driven by this, at UT22:00, in the case of the electron density reduction in other inversion area, the electron density at 350 km from 45N to 55N increases comparing with the UT18:00 image. This distinctive increase phenomenon proves the feasibility of forecasting the 3-D ionosphere change with 3-D velocity images. We will verify this phenomenon with ionosonde data next. At UT22:00, the electron density peak goes back to the 350-km height, from Fig. 8(l), at 250–350 km, the electron density decreases in maximum strength of the day, but the changes of other ionosphere layers are relatively smooth. On the whole, taking 4:00 hour’s time resolution into account, the electron density variation is in good agreement with the velocity images, and it shows that this study presents a
way to forecast the 3-D electron density variation with the 3-D velocity images. Furthermore, we use the conventional algorithm MART to reconstruct the ionosphere electron density images at the same area and the same time and compare the MART’s results with the new algorithm’s results among different latitude surfaces. The comparisons are shown in Table I. In the table, it is shown that with the latitudes reduced, the differences of maximum value, averages, and standard deviations are gradually increasing; specifically, the maximum value is increased by one order of magnitude. However, considering Fig. 8(g), the electron density peak value reaches 1012 el/m3 at 14:00UT, while the differences are much less than the magnitude. That means that the inversion results of the two algorithms are basically consistent. In addition, the ionosphere ionosonde data are downloaded and used to prove the new algorithm’s results. As aforementioned, the new algorithm detects that, on 12◦ E surface at UT18:00, between 45◦ N and 55◦ N, there is an electron density increase phenomenon at 300–400 km, which is inconsistent with the overall change. This paper selects ionosonde data at European regions to verify the reliability of this phenomenon. In this area, there are three ionosonde observatory stations, but only JR055 JULIUSRUH (54.60◦ N, 13.40◦ E) conducted ionosphere exploration that day. In Fig. 8(j), the electron density increases at 300–400 km and decreases at 100–200 km, respectively. JR055 conducts 11 measurements from UT17:28 to UT20:13. We intercept 100- to 200-km and 300- to 400-km detection results in Fig. 9. In the figure, the IRI model’s results and the new algorithm’s results (the UpIRI) are also printed out as the reference.
2382
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 53, NO. 5, MAY 2015
Fig. 10. Electron density profiles of different models at 54.60◦ N, 13.40◦ E. The top panel is the electron density profile at UT14:00, whereas the bottom panel is at UT20:00.
As shown in the figure, the detection results of ionosonde have exactly the same variation process comparing with the new algorithm’s results. Particularly on layers 300–400 km, the IRI model fails to reflect the real electron density increasing phenomena, whereas the new algorithm properly reflects the abnormal disturbance with high electron density detection accuracy. The velocity images of the new algorithm on UT18:00 are consistent with the ionosonde’s detection results. Fig. 10 shows the electron density profiles of different models at the JR055 station site. In the top panel in Fig. 10, the electron density peak height experiences a decline, comparing the ionosonde data with the IRI model (from 350 to 300 km). The new algorithm result shows the peak height decline exactly. The bottom panel gives the comparison at UT20:00. The electron density value of ionosonde data is greatly increased comparing with the IRI model on layers 200–500 km, whereas the new algorithm gives more precise results. V. C ONCLUSION The real-time detection and analysis of the Earth’s space environment abnormal outbreak activity caused by some natural disasters are likely to be a potential and ancillary means to forecast major natural disasters. The ionosphere tomography technology based on GNSS can obtain the 3-D ionosphere morphology. The inversion algorithm is the critical challenge, which restricts the precision and accuracy of the ionosphere tomography technology. In this paper, by adding the loose constraints among the grids,
forming a stable mesh structure and then composing the virtual normal equation, combining the virtual normal equation with the observed normal equation, the electron density parameters are estimated with best fit to the real observations based on the priori information. In addition, different from previous algorithms, the new algorithm introduces the electron density velocity parameters, which is more reasonable with actual ionosphere variation. The 3-D velocity profiles are useful and provide multisource information for studying the spatial ionosphere variation, particularly the ionosphere variation under ionosphere disturbances. In addition, the 3-D ionosphere changes velocity can provide reliable ionosphere gradient information in each layer for some joint-inversion algorithms and then improve the accuracy of ionosphere model. We have used the new algorithm to reconstruct the ionosphere electron density and velocity images on August 15, 2003, around Europe. By comparing the results with the conventional algorithm MART’s results and ionosphere ionosonde data, the feasibility and reliability of the new algorithm are verified. The related research provides a new way to real-time detection and prediction of ionosphere change. ACKNOWLEDGMENT The authors would like to thank IGS for providing GNSS data and global ionospheric radio observatory (GIRO) for providing the ionosonde data. The authors would also like to thank the anonymous reviewers for their constructive comments in improving this paper.
YAO et al.: IONOSPHERE TOMOGRAPHY ALGORITHM WITH CONSTRAINTS AND VELOCITY PROFILE
R EFERENCES [1] M. M. Alizadeh, H. Schuh, S. Todorova, and M. Schmidt, “Global ionosphere maps of VTEC from GNSS, satellite altimetry, formosat3/COSMIC data,” J. Geodesy, vol. 85, no. 12, pp. 975–987, Dec. 2011. [2] C. K. Avinash and S. Malcolm, Principles of Computerized Tomographic Imaging. Piscataway, NJ, USA: IEEE Press, 1988, pp. 275–296. [3] K. Davies and D. M. Baker, “Ionospheric effect observed around the time of the Alaskan earthquake of March 28, 1964,” J. Geophys. Res., vol. 70, no. 9, pp. 2251–2253, May 1965. [4] R. S. Leonard and R. A. Barnes, “Observation of ionospheric disturbances following the Alaska earthquake,” J. Geophys. Res., vol. 70, no. 5, pp. 1250–1253, Mar. 1965. [5] D. Dettmering, M. Schmidt, and R. Heinkelmann, “Combination of different space-geodetic observations for regional ionosphere modeling,” J. Geodesy, vol. 85, no. 12, pp. 989–998, Dec. 2011. [6] R. Garcia and F. Crespon, “Radio tomography of the ionosphere: Analysis of an underdetermined, ill-posed inverse problem, regional application,” Radio Sci., vol. 43, no. 2, pp. RS2014-1–RS2014-13, Apr. 2008. [7] T. Hobiger, T. Kondo, and Y. Koyama, “Constrained simultaneous algebraic reconstruction technique (C-SART)—A new and simple algorithm applied to ionospheric tomography,” Earth Planets Space, vol. 60, no. 7, pp. 727–735, Jul. 2008. [8] Z. Z. Liu, “Ionosphere tomographic modeling and applications using Global Positioning System (GPS) measurements,” Univ. Calgary, Calgary, AB, USA, 2004. [9] J. Y. Liu et al., “Pre-earthquake ionospheric anomalies registered by continuous GPS TEC measurements,” Ann. Geophys.-Germany, vol. 22, no. 5, pp. 1585–1593, 2004. [10] E. P. Macalalad, L. C. Tsai, J. Wu, and C. H. Liu, “Application of the Taiwan Ionospheric Model to single-frequency ionospheric delay corrections for GPS positioning,” GPS Solutions, vol. 17, no. 3, pp. 337–346, Jul. 2013. [11] S. Kon, M. Nishihashi, and K. Hattori, “Ionospheric anomalies possibly associated with M ≥ 6.0 earthquakes in the Japan area during 1998–2010: Case studies and statistical study,” J. Asian Earth Sci., vol. 41, no. 4/5, pp. 410–420, Jun. 2011. [12] A. Krankowski, I. Zakharenkova, A. Krypiak-Gregorczyk, I. I. Shagimuratov, and P. Wielgosz, “Ionospheric electron density observed by FORMOSAT-3/COSMIC over the European region and validated by ionosonde data,” J. Geodesy, vol. 85, no. 12, pp. 949–964, Dec. 2011. [13] P. F. Fougere, “Ionospheric radio tomography using maximum entropy,” Radio Sci., vol. 30, no. 2, pp. 429–444, Mar./Apr. 1995. [14] S. E. Pryse and L. Kersley, “A preliminary experimental test of ionospheric tomography,” J. Atmos. Terrestrial Phys., vol. 54, no. 7/8, pp. 1007–1012, Jul./Aug. 1992. [15] S. A. Pulinets and K. Boyarchuk, Ionospheric Precursors of Earthquakes. Berlin, Germany: Springer-Verlag, 2004. [16] M. Schmidt, D. Bilitza, C. K. Shum, and C. Zeilhofer, “Regional 4D modeling of the ionospheric electron density,” Adv. Space Res., vol. 42, pp. 782–790, 2008. [17] S. Saroso, J. Y. Liu, K. Hattori, and C. H. Chen, “Ionospheric GPS TEC anomalies and M > 5.9 earthquakes in Indonesia during 1993–2002,” Terrestrial Atmos. Ocean. Sci., vol. 19, no. 5, pp. 481–488, Oct. 2008. [18] S. Todorova, H. Schuh, and T. Hobiger, “Using the global navigation satellite systems and satellite altimetry for combined global ionosphere maps,” Adv. Space Res., vol. 42, no. 4, pp. 727–736, Aug. 2007. [19] L. C. Tsai, W. H. Tsai, W. S. Schreiner, F. T. Berkey, and J. Y. Liu, “Comparisons of GPS/MET retrieved ionospheric electron density and ground based ionosonde data,” Earth Planets Space, vol. 53, pp. 193–205, 2001.
2383
[20] D. B. Wen and Y. B. Yuan, “A hybrid reconstruction algorithm for 3-D ionospheric tomography,” IEEE Trans.Geosci. Remote Sens., vol. 46, no. 6, pp. 1733–1739, Jun. 2008. [21] D. B. Wen and Y. Wang, “A new two-step algorithm for ionospheric tomography solution,” GPS Solutions, vol. 16, no. 1, pp. 89–94, Jan. 2012. [22] Y. B. Yao et al., “Analysis of pre-earthquake ionospheric anomalies before the global M = 7.0+ earthquakes in 2010,” Nat. Hazards Earth Syst. Sci., vol. 12, no. 3, pp. 575–585, 2012. [23] X. N. Yue, W. S. Schreiner, and Y. H. Kuo, “Evaluating the effect of the global ionospheric map on aiding retrieval of radio occultation electron density profiles,” GPS Solutions, vol. 17, no. 3, pp. 327–335, Jul. 2013.
Yibin Yao received the B.Sc., Master’s, and Ph.D. degrees (with distinction) in geodesy and surveying engineering from Wuhan University, Wuhan, China, in 1997, 2000, and 2004, respectively. He is currently a Professor with Wuhan University. His main research interests include global navigation satellite system ionospheric/atmospheric/ meteorological studies, theory and method of surveying data processing, and GPS/MET and highprecision GPS data processing. He has authored over 50 publications since 2001 and six books in those areas. Prof. Yao has been a recipient of a Luojia Professorship at Wuhan University, supported by the Program for New Century Excellent Talents in University.
Jian Kong received the B.Sc. degree from Shandong University of Science and Technology, Qingdao, China, in 2009. He is currently working toward the Ph.D. degree in geodesy and surveying engineering in the School of Geodesy and Geomatics, Wuhan University, Wuhan, China. His main research interests include GNSS data processing, global ionospheric modeling using multisource geodesy observations, 3-D ionospheric tomography based on GNSS observations, and ionospheric anomalies analysis under abnormal conditions.
Jun Tang received the B.Sc. degree in surveying and mapping engineering and the M.Sc. degree in geodesy and surveying engineering from East China Institute of Technology, Fuzhou, China, in 2006 and 2009, respectively. He is currently working toward the Ph.D. degree in geodesy and surveying engineering in the School of Geodesy and Geomatics, Wuhan University, Wuhan, China. His research interests include ionospheric modeling using diverse data, multidimensional ionospheric tomography based on global navigation satellite system (GNSS) observations, and GNSS applications.
