Assessment of ionosphere tomographic modeling performance ... - Wiley

RADIO SCIENCE, VOL. 41, RS1007, doi:10.1029/2004RS003236, 2006

Assessment of ionosphere tomographic modeling performance using GPS data during the October 2003 geomagnetic storm event Zhizhao Liu,1,2 Susan Skone,1 and Yang Gao1 Received 15 December 2004; revised 3 October 2005; accepted 10 November 2005; published 22 February 2006.

[1] Precise ionospheric modeling is important for single-frequency Global Positioning

System (GPS) users to achieve optimal positioning accuracy because the ionospheric signal delay is now the largest error source for positioning and navigation with GPS. The ionospheric modeling during ionospheric storms is particularly critical since the signal delay may be higher than normal and may differ significantly from the broadcast ionosphere model (currently employed by single-frequency users). In this study, a tomographic technique is used to model the ionosphere over North America using data collected from a network of dual-frequency GPS receivers. In support of real-time applications of the ionosphere model, short-term (5-min) forecasts of ionospheric total electron content (TEC) are also performed. To validate the accuracy of the forecast ionospheric TEC, a comparison of the TEC predictions with the observed TEC data (which are inferred from dual-frequency GPS observations) is carried out. Analyses are conducted using GPS data recorded during a 2003 geomagnetic storm event (29– 31 October). Results indicate that under less disturbed conditions, an average accuracy of 5 6.5 total electron content units (TECU, 1 TECU = 1016 el m2) can be obtained for the vertical TEC prediction and that 80% of slant TEC can be recovered by the model predictions. During extreme ionospheric storm periods (Kp = 9), the vertical TEC forecasting accuracy has a degradation of 2 3 TECU from the 3-day mean value, and the relative error is several percent to 10% larger than the 3-day average level. Citation: Liu, Z., S. Skone, and Y. Gao (2006), Assessment of ionosphere tomographic modeling performance using GPS data during the October 2003 geomagnetic storm event, Radio Sci., 41, RS1007, doi:10.1029/2004RS003236.

1. Introduction [2] The global positioning system (GPS) signals on two frequencies (L1 = 1575.42 MHz and L2 = 1227.60 MHz) provide a convenient way to precisely determine the amount of total electron content (TEC) along the GPS signal path, given that range observations on both frequencies are recorded by GPS receivers. Therefore the refraction delay in GPS signals introduced by ionosphere could be readily compensated. However, this approach is not feasible for single-frequency

1 Department of Geomatics Engineering, University of Calgary, Calgary, Alberta, Canada. 2 Now at Leica Geosystems, Ltd., Calgary, Alberta, Canada.

Copyright 2006 by the American Geophysical Union. 0048-6604/06/2004RS003236

receivers where the TEC cannot be estimated directly using the dispersive nature of the ionosphere. An ionospheric modeling technique is hence required for all GPS users with single-frequency receivers. The employment of ionospheric models allows single-frequency users to correct the ionospheric refraction error and consequently to improve their positioning and navigation precision. Because of the significant uncertainty of ionospheric electron distribution in both space and time, precise modeling of ionospheric refraction delay is still a challenging task. A broadcast ionosphere model is currently available for single-frequency users [Klobuchar, 1986], as transmitted in the GPS navigation message. This model is simple, however, being based on eight coefficients to describe the spatial distribution of TEC, it mitigates only 50% of the overall ionospheric effect [Klobuchar, 1987]. An alternative technique is to model the ionosphere on the basis of dual-frequency observations from a permanent network of GPS receivers.

RS1007

1 of 12

RS1007

LIU ET AL.: IONOSPHERE MODELING USING 2003 GPS DATA

Observations of slant TEC are derived for all satellitereceiver lines of sight within the network, and a spatial model of the vertical TEC or three-dimensional electron density is derived for the given region. This approach is employed and tested here. [3] Ionospheric modeling is particularly difficult during periods of increased solar activity, where associated disturbances may result in large gradients in TEC [e.g., Coster et al., 2003]. The October – November 2003 solar event studied in this paper was initiated by three very large sunspot clusters [Doherty et al., 2004]. The sudden increase of sunspots resulted in a series of solar flares and coronal mass ejections [Department of Commerce, 2004]. These events caused a severe geomagnetic storm in the Earth’s environment and disturbance in the ionosphere [Doherty et al., 2004]. The ionospheric disturbance was directly reflected in the large TEC variations, which in turn led to significant range delay or advance in the GPS signals. Extraordinary solar storms were observed during 18 October to 5 November 2003, with 44 M-class flares and 11 X-class flares [Woods et al., 2004]. These significant solar storm events are detrimental to GPS users. Ionospheric modeling with GPS measurements may be difficult because of reduced data quality, particularly for the L2 signal [Skone, 2001]. The loss of L2 signal will disable the derivation of carrier phase TEC data from L1 and L2 carrier phase measurements. Thus the TEC smoothing procedure in the ionospheric modeling which uses the phase-derived TEC data to smooth code-derived TEC cannot be performed because of the absence of L2 signal at a given epoch. At those particular epochs, when L2 signal is lost, however, code-derived TEC data can be used alone in the ionospheric modeling because pseudorange measurements are always available. Thus the effect of losing L2 signal on modeling is not significant though the code-derived TEC data have a lower accuracy than the carrier phase – smoothed one. Without a sufficiently dense reference network, or adequate modeling techniques, localized TEC gradients may not be resolved, and ionosphere model accuracies may be degraded. Such storms provide excellent opportunities, however, for testing of ionospheric modeling methods under challenging conditions. In this investigation, the performance of a tomographic ionospheric model was analyzed using GPS data recorded at a high level of ionospheric activity (Kp = 8 9) during 29– 31 October 2003. In these 3 days, the Kp indices reached maximum values of 9 over a total period of 18 hours.

2. Model Description [4] The ionospheric model implemented in this investigation is the tomographic model [Raymund et al., 1990;

RS1007

Hansen et al., 1997; Howe et al., 1998; Liu and Gao, 2003]. Generally, there are two types of approaches that could be categorized as tomographic models, namely, the voxel-based approach and the function-based approach [Liu, 2004]. By employing the function-based approach, the electron distribution in the ionosphere is described by two functions, namely, spherical harmonic functions for the horizontal property and empirical orthogonal functions (EOF) for vertical property. The voxel-based approach assumes that the ionosphere comprises numerous small cells, each of which has a uniform electron density. In the voxel-based approach, the cell size must be chosen appropriately before modeling to ensure that this assumption is as close to reality as possible. In practical modeling implementations, the number of voxels usually is very large even in a small area. For example, a 4 4 40 voxel grid was defined over a region of 400 km2 and 15 km in height to model the troposphere [Flores et al., 2000]. To use the voxel-based approach to model ionospheric electron density, more voxel cells are required in order to achieve the same resolution because the ionosphere is much thicker than the troposphere. In such a way, intensive computations will be required to perform the estimation of unknown parameters, and it will be computationally difficult in real-time implementation. In comparison, the function-based approach has significantly less unknown parameters in the modeling. The determination of total unknown ionospheric parameters is dependent on the number of vertical ionospheric layers (denoted by K) and highest order of spherical harmonic functions (denoted by M). The total number of parameters can be calculated as K(M + 1)(2M + 1) [Liu and Gao, 2004]. The selection of the number of unknown parameters is a function of several factors, including the size of ionospheric modeling area, the expected modeling accuracy, computational burden, and the density of GPS data. When the GPS satellite and receiver interfrequency biases are estimated, the number of parameters will increase by (nr + ns), where nr and ns represent the number of GPS receivers and GPS satellites tracked by the GPS receivers, respectively. Our data analysis shown below indicates that the functionbased approach requires a much smaller number of parameters in the ionospheric modeling than the voxelbased approach. Thus the function-based tomographic approach is chosen for this investigation. Assuming that the a priori ionospheric electron density function N0(l, f, z) may be obtained from other ionospheric models or empirical ionospheric data and that the correction part to N0(l, f, z) is modeled by a function-based tomographic approach, the mathematical model can be described as follows [Liu and Gao, 2003]. N0(l, f, z) is just an approximate value of Ne(l, f, z), and it can be obtained by many methods, for instance, from a parameterized ionospheric model, or derived from data of GPS occul-

2 of 12


RS1007

Bp satellite’s differential instrumental delays on code pseudorange measurements between the L1 and L2 frequencies; g s qu ar e d L 1 a n d L2 f re q ue nc y r a t i o ; f1 2 1575:42 2 77 2 ¼ ¼ . g¼ 1227:6 60 f2

tation, ionosonde, or incoherent scatter radar [Daniell et al., 1995; Hajj and Romans, 1998]:

dNe ðl; f; zÞ ¼ þ

bm nk

K M M X X X am nk cosðmlÞ

k¼1 m¼M n¼jmj m sinðmlÞ Pn ðcos fÞZk ðzÞ

ð1Þ

where dNe(l, f, z) correction to the a priori electron density function N0(l, f, z); (l, f, z) three coordinate components of the spatial position, representing Sun-fixed longitude, geomagnetic latitude, and altitude, respectively; Pm n (cos f) associated Legendre polynomial of order m and degree n(0 m n); Zk(z) empirical orthogonal functions (EOF); model’s coefficients to be estimated; am nk model’s coefficients to be estimated; bm nk K highest order of empirical orthogonal functions; M highest order of spherical harmonics functions. The summation of correction part dNe(l, f, z) and the a priori electron density function N0(l, f, z) forms a good estimation of the ionospheric electron density function, denoted by Ne(l, f, z). With knowledge of Ne(l, f, z), the ionospheric total electron content along any line of sight can be derived through integration. These modeled TEC values can then be used by single-frequency GPS users to correct for the ionospheric effect. [5] In the ionospheric modeling, TEC observations are first obtained from dual-frequency GPS data in a reference network, and dNe(l, f, z) is then estimated using the TEC measurements. For instance, the TEC data can be derived as below when adopting code pseudorange GPS measurements:

TECR ¼

f 21 ½ðP1

p

P2 Þ B i B 40:3ð1 gÞ

ð2Þ

The TECR data derived from code pseudorange GPS measurements is quite noisy because of a high level of noise in pseudorange measurements P 1 and P2. A smoothing algorithm [Skone, 1998; Liu, 2004] with the aid of much more precise carrier phase GPS measurements on both L1 and L2 frequencies is employed in this investigation to reduce the noise in the TECR measurements. The carrier phase – smoothed TEC data, denoted by TECSM, has a relationship with dNe(l, f, z) and N0(l, f, z) that can be expressed as

TECSM ¼

¼

Zsv rx Zsv

½N0 ðl; f; zÞ þ dNe ðl; f; zÞds

N0 ðl; f; zÞds þ

rx

Zsv

dNe ðl; f; zÞds

rx

¼ TEC0 þ

Zsv

dNe ðl; f; zÞds

ð3Þ

rx

After the TEC0 term, which is derived from the a priori electron density function N0(l, f, z), is subtracted from both sides of equation (3), the following expression is derived: dTEC ¼ TECSM TEC0 ¼

Zsv

dNe ðl; f; zÞds

ð4Þ

rx

Inserting equation (1) into equation (4),

dTEC ¼

K M M X X X

am nk

k¼1 m¼M n¼jmj

þ

where

K X

M X

M X

Zsv

m

m

cosðmlÞPn ðcos fÞZk ðzÞds

rx

bm nk

k¼1 m¼M n¼jmj

TECR TEC measurement derived from GPS code pseudorange; P1 code pseudorange measurement on L1 frequency; P2 code pseudorange measurement on L2 frequency; Bi receiver’s differential instrumental delays on code pseudorange measurements between the L1 and L2 frequencies;

RS1007

Pn ðcos fÞZk ðzÞds

Zsv

sinðmlÞ

rx

ð5Þ

Equation (5) is the fundamental analytical expression of the function-based tomographic approach. The unknown m parameters in this approach are coefficients am nk and bnk. p Interfrequency biases Bi and B are also implicitly

3 of 12

RS1007


included in dTEC. In the parameter estimation, these two types of bias parameter should be moved to the righthand side of equation (5), and they are estimated along m with the coefficients am nk and bnk. The estimator used in this study is a Kalman filter, which allows real-time processing of GPS measurements and performing tomographic model construction in real time [Liu, 2004]. The Kalman filter estimator can be described by the following two equations [Brown, 1983]: xk ¼ Fk;k1 xk1 þ wk1

ð6aÞ

z k ¼ H k xk þ vk

ð6bÞ

where xk system state vector consisting of the unknown parameters (model coefficients and interfrequency biases) at time tk; Fk,k1 transition matrix relating the estimated xk1 and predicted xk; wk1 system noise vector assumed to be white uncorrelated sequence with known covariance; zk measurement vector at time tk; Hk design matrix describing the relationship between measurement and state vector at time tk; vk measurement error vector assumed to be a white sequence with known covariance. The equation (6a) is the dynamic equation, and the (6b) is called the measurement (observation) equation [Brown, 1983]. The equation (6a) describes a predictive relationship between the current system state and the future state. In the dynamic equation, the parameters in vector xk, namely, ionosphere model coefficients, are modeled as a first-order Gauss-Markov process. The corresponding transition matrix in the system model can be written as 2 b Dt 3 0 ... 0 e 1 6 0 eb2 Dt . . . 0 7 7 ð7Þ Fk;k1 ¼ 6 4 ... ... ... ... 5 0 0 . . . ebn Dt and the process noise matrix Qk can be written as

s21 1 e2b1 Dt

0 2b Dt 2 6 2 0 s 2 1e Qk ¼ 6 4 0 0 2

3 0 7 0 7 5

2 2bn Dt sn 1 e ð8Þ

In equations (7) and (8), s and 1/b are standard deviation of random process noise and correlation time of the process, respectively. During ionospheric-disturbed per-

RS1007

iods, the ionospheric TEC variation will be significant compared to the quiet periods, but the variation of the ionosphere can be properly modeled by an appropriate specification of the standard deviation and correlation time in equation (8). When the dynamics of the ionosphere variation increases, the standard deviation should be increased to account for the increasing disturbance in the dynamic model, and the correlation time should be accordingly shortened. In the ionospheric TEC prediction, the predictive functionality of the dynamic equation will be employed.

3. Analysis [6] As stated above, this investigation focuses on the study of ionospheric model performance for a functionbased tomographic approach under a high level of ionospheric activities (Kp = 8 9). For such purposes, 3 days of GPS measurements recorded during 29– 31 October 2003 are analyzed. During the 3 days, extraordinary geomagnetic activities were observed, and a total of 18 hours with Kp index value of 9 were experienced within that period [Woods et al., 2004; Komjathy et al., 2004]. The Kp values over the 3 days are depicted in Figure 1. The analysis of GPS measurements during that period therefore will sufficiently serve the purpose of testing the performance of the function-based tomographic model under some of the most challenging conditions. This study focuses on storm effects in North America; strong gradients during periods of stormenhanced density were observed across North America during this event [Komjathy et al., 2004]. During geomagnetic storms, a phenomenon known as stormenhanced density, which is driven by processes in Earth’s magnetosphere and is associated with large gradients in the ionospheric and plasmaspheric TEC, is often observed [Coster et al., 2003]. Large gradients in TEC are primarily observed in the equatorial and polar regions, but it can also be observed at midlatitudes during moderate to severe geomagnetic disturbances [Coster et al., 2003]. In the October 2003 event, stormenhanced densities were observed over the continental United States during 2000 UT on 29 October to 0100 UT on 30 October and were observed again from 1800 UT on 30 October to 0200 UT on 31 October [Doherty et al., 2004]. The GPS measurements were collected from 55 International GNSS Service (IGS) stations distributed in North America. The geographical distribution of the GPS stations is shown in Figure 2. The coverage of the stations is from 25.73N to 82.49N in latitude and 52.68W to 147.50W in longitude. Among the 55 stations, only GPS measurements of 50 stations are used to construct the tomographic model, and the remaining 5 stations are used as independent ‘‘truth’’ data stations to calibrate the modeling accuracies. These five stations are

4 of 12

RS1007


RS1007

into the Kalman filter estimator. The interfrequency biases estimated from the first run are applied as known m values, and only the model coefficients am nk and bnk are m estimated. Once the ionospheric coefficients ank and bm nk are estimated, they can be broadcast to all model users, such as single-frequency GPS users, to calculate TEC values. 3.2. Model Accuracies

Figure 1. Kp values during 29– 31 October 2003. AMC2, ERLA, HRN1, IDTD, and PCK1, and they are circled in Figure 2. All 55 stations are equipped with dual-frequency GPS receivers and can output P1 and P2 code pseudorange as well as F1 and F2 carrier phase measurements on L1 and L2 frequencies. This allows the derivation of smoothed TEC data from GPS measurements, and the carrier phase– smoothed TEC data serve as input data for the ionosphere tomographic model. All 55 IGS stations record GPS measurements at a rate of 30 s. This rate is also adopted in the ionospheric modeling. The elevation cutoff angle used for editing GPS measurement during both modeling and calibration is 15. 3.1. Model Construction [7] Before the accuracy estimation is performed at the 5 test stations, the ionospheric model has to be constructed for the area covered by the 50 stations. In this study, 15 min of GPS measurements from the 50 stations are used to construct the tomographic model as described in the previous section. TEC data are derived from dualfrequency GPS measurements collected at each epoch during the 15-min period, and they are then input into a Kalman filter. The Kalman filter estimates the ionom spheric coefficients am nk and bnk as well as receiver and satellite interfrequency biases Bi and Bp for each GPS receiver and satellite. During the modeling, the software is run twice. In the first run, GPS measurements from all 55 stations are input into the Kalman filter estimator. The purpose of the first run is to estimate the interfrequency biases Bi and Bp for all 55 GPS stations and the tracked satellites. The interfrequency estimates are then used as known values in both modeling and calibration procedures. In the second run, only the measurements from the 50 stations designated for model construction are input

[8] According to equation (5), the correction part of the TEC value, dTEC, at any user station can be calculated by using the ionospheric coefficients once the broadcast coefficients are received. The summation of dTEC and the a priori value TEC0 forms the total electron content value at a given user station. The TEC value is an important quantity that is used in determining the magnitude of ionospheric delay on GPS signals. One crucial element that represents the performance of an ionospheric model is how well the calculated TEC data agree with the truth TEC data, which can be obtained from independent observations or other models. In this study, the calibration truth TEC data are derived from dual-frequency GPS measurements. As stated previously, the modeling calibration is performed at five test stations. The 5 test stations are assumed to be independent of the other 50 stations, considering that their GPS measurements have not been used in the modeling construction. However, it should be noted that the GPS receiver interfrequency biases for these five calibration stations are obtained from the first run in estimating all the GPS receiver and satellite interfrequency biases together. Thus a certain correlation between the interfrequency biases for the five stations and the rest of the stations may exist during the first run. Considering that the magnitude of interfrequency biases are much smaller than

Figure 2. Geographical distribution of GPS stations.

5 of 12


RS1007

that of the ionospheric TEC values and to simplify the TEC calibration procedure, this type of correlation is considered ignorable during the TEC calibration. Taking advantage of the dual-frequency GPS measurements, smoothed TEC observations are computed for the five test stations, and these TEC data are used as a reference with which the calculated TEC predictions are compared. [9] Two indicators are used in this study to assess the agreement of the calculated TEC data with the truth calibration TEC data. One is a relative indicator which denotes the disagreement in percentage of the calculated TEC with respect to the truth TEC data. The other one is an absolute indicator which statistically calculates the root-mean-square (RMS) of the differences between the calculated TEC and truth TEC data after mapping them to the zenith direction. The calculated TEC data and truth TEC values are denoted as dc and do, respectively. It is worth noting that in the space domain, the calculation of TEC data is performed at user stations, namely, five test stations in this study and that in the time domain, the calculation is performed for the future 5-min epoch at the user stations. Therefore the calculation in fact represents a short-term prediction in time and an interpolation in space. The realization of a short-term (5-min) prediction is crucial in practice for real-time TEC data users, to allow buffer time for calculation, transmission, and implementation, such as GPS real-time kinematic positioning and navigation, or the Wide Area Augmentation System (WAAS) [Liu and Gao, 2003]. The difference between calculated and measured TEC data is given as Dd ¼ dc do

ð9Þ

where Dd is considered to be the model error. The magnitude of this difference shows the agreement level of the calculated TEC data with respect to the calibration (observed) TEC. Small differences reflect a good agreement between the two sets of TEC values. If the differences are further mapped to the zenith direction, the vertical TEC (VTEC) errors are obtained. One VTEC error can be calculated for each comparison. For all the TEC comparison results, a root-mean-square (RMS) value can be estimated on the basis of the VTEC errors. Thus the absolute accuracy indicator, RMSvtec, can be defined as

RMSvtec ¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uP un u Ddvtec;i Ddvtec;i t i n

RS1007

zenith direction has an advantage of being able to eliminate the dependence of the TEC values on the elevation angles, which is more convenient for the consistent evaluation of ionosphere model accuracies. In contrast, the relative indicator, named relative error (RE), can be defined as RE ¼

n 1X Ddi =do 100% i n i¼1

ð11Þ

where Ddi is the ith modeling error, doi is the ith TEC observation, and n denotes the total number of TEC comparisons. The relative error represents the mean percentage of ionospheric modeling error that the calculated TEC data may contain with respect to the observed reference ionospheric TEC. For a good ionospheric model, this relative error should be as small as possible. In comparison, the VTEC RMS indicator can be treated as an ‘‘absolute’’ indicator to describe the accuracy of the calculated TEC data, while the relative error is a ‘‘relative’’ indicator. It can be observed that the computation of relative error is affected by the magnitude of the observation TEC data. Given a certain Ddi, the relative error is determined by the size of doi ; however, this is not the case for VTEC RMS estimation. [10] Since the calibration TEC observations are directly derived from dual-frequency GPS measurements at the five test stations, they inherently contain the satellite and receiver interfrequency biases. Nevertheless, the predicted TEC data do not include the interfrequency biases because they are computed from the tomographic model. Thus it is necessary to cancel out the effect of the interfrequency biases in the calibration TEC when calculating in the difference Dd. It should be noted that the estimated values of the receiver and satellite interfrequency biases have been obtained from the first run of the modeling software. Thus the estimated biases can be added to predicted TEC data. Therefore the difference Dd should be free of bias effect, and it contains only the modeling error with respect to the truth TEC data. In the modeling calibration procedure, differences between predicted and calibration TEC data are examined. In section 4, both VTEC RMS error and relative error for the calculated TEC data at five test stations will be presented.

4. Results ð10Þ

where Ddvtec,i is the value of the ith modeling error after mapping to the zenith direction. The projection to the

[11] GPS receiver tracking performance under geomagnetic storms may degrade significantly compared to normal geomagnetic conditions [Skone, 2001]. In the data compilation for deriving smoothed TEC data, it is found that GPS carrier phase measurements suffer from

6 of 12

RS1007


RS1007

Figure 3. Comparison between the predicted and truth TEC data. Blue pluses represent the truth TEC data observed at the test stations and red dots denote the predicted TEC data for the test stations.

severe cycle slips. Cycle slip detection and repairing techniques [Blewitt, 1990; Bisnath, 2000] are employed to correct the GPS carrier phase observations; however, there are observation epochs when the correct repairing of cycle slips were not possible. At such observation epochs, the carrier phase – smoothing algorithm is not implemented since carrier phase measurements at those epochs are not usable because of irreparable cycle slip or carrier phase data loss. However, pseudorange code measurements are always available and free of cycle slip. In such cases, the code-based TEC data are used as TEC measurements for modeling. When the code-based TEC data are used, the change of observation accuracy has to be accordingly considered in the modeling since the code-based TEC data have a lower accuracy than the smoothed one. In the model construction, only 15 min of GPS data are used for each session. Thus this approach is in effect for a short period (15 min), and the degradation of TEC measurement noise is limited. Using this approach, the continuity of TEC data is ensured. In each new session, a new offset between carrier phase and pseudorange TEC data is calculated,

and the cycle slip detection and repairing procedure commences again. [12] The modeling was performed using 15-min GPS measurements from the 50 GPS stations, and the 5-min temporal prediction was performed at the 5 test stations. The predicted TEC data were compared to the truth TEC data that were directly derived from GPS measurements of the five test stations. A sample of the TEC comparison is shown in Figure 3, where blue pluses represent the truth TEC data observed at the test stations and red dots denote the predicted TEC data for the test stations. At each test station, only the satellite with the maximum amount of tracking time was selected for the comparison. Figure 3 shows that the TEC data predicted by the constructed tomographic model for each test station had a good agreement with the truth data observed at the same test station. [13] The modeling days are 29 – 31 October 2003, during which there are two periods (period I, 1900 UT, 29 October to 0300 UT, 30 October; and period II, 1900 –2400 UT, 30 October) have high value of Kp index (Kp = 9) as shown in Figure 1. Figure 4 shows that

7 of 12

RS1007


RS1007

Figure 4. TEC variation observed at five test stations during 29– 31 October 2003. significant TEC increases were observed at all the test stations during both period I and period II. Because of incomplete data at station AMC2 on 31 October 2003, the AMC2 station is excluded from analysis for that day. Since the major objectives of this study are to investigate the ionospheric prediction accuracy during storm periods and the impact of storm enhanced density on the prediction accuracy, the TEC comparison statistical results obtained during the storm periods are presented in Figures 5 and 6. Figure 5 shows the vertical RMS error on the basis of hourly statistical results. In each hour, there were 12 sessions of 5-min TEC prediction data. The prediction data were compared to the truth TEC data, and the vertical RMS error could be calculated from the TEC differences as equation (10) showed. The vertical RMS error was calculated on the basis of every 12 sessions of TEC comparison results. Figures 5 and 6 show that the hourly prediction results had a good accuracy. For the vertical RMS errors, they generally were varying between 4.0 6.0 total electron content units (TECU, 1 TECU = 1016 el m2), and for the relative errors, they were varying between 20% 25% during most periods at the test stations, except for the two periods (period I, 1900 UT, 29 October to 0300 UT,

30 October; and period II, 1900 –2400 UT, 30 October). During the two periods, the vertical RMS errors were much larger. Figure 5 showed that the hourly prediction errors were in the range 6 10 TECU during periods I and II, which was significantly higher than during other periods. The maximum prediction error reaches 11.3 TECU at AMC2 station, as shown in the top plot of Figure 5. Examining Figure 1, it is found that the higher Kp indices were observed during the two periods. Two major storm events were observed during these two periods [Komjathy et al., 2004]. Higher values of Kp index usually indicate that ionospheric activity is more dynamic and that it is more difficult to resolve local ionosphere features using a model. Similar modeling results were also obtained by using a WAAS planar fit approach [Komjathy et al., 2004]. Results showed the RMS error for the slant planar fit residuals was 2.6 m for 29 October and 3.4 m for 30 October 2004, which were approximately equivalent to 16.0 and 20.9 TECU in the slant, respectively [Komjathy et al., 2004]. Assuming that the average elevation angle of the TEC measurements in planar fitting was 30, the factor of converting slant residual to the vertical was about 0.588 when using a single-layer mapping function. After mapping the slant

8 of 12

RS1007


RS1007

Figure 5. Vertical TEC prediction RMS error at test stations.

residual to the vertical, the residual RMS errors were 9.4 and 12.3 TECU for 29 and 30 October 2004, respectively. It should be noted that the vertical TEC RMS errors studied in this investigation were calculated from a comparison with external observation truths. [14] The average TEC prediction errors for the 3-day period as well as period I and period II are summarized in Table 1. Compared to the 3-day average TEC prediction accuracy, a significant deterioration of the prediction accuracy could be seen for both period I and period II. During the two periods, the VTEC predictions had larger errors than the 3-day mean VTEC prediction errors by 2.0 3.0 TECU, as shown in Table 1. This was equivalent to an accuracy degradation of 38% for period I (averaging among the five test stations) to 48% for period II (averaging among the five test stations) from the 3-day mean accuracy. [15] Figure 6 shows the relative errors calculated using prediction and truth TEC data for each station. The

calculation of relative errors was also based on TEC comparison data for each 1-hour period. The mean relative errors are summarized in Table 1, which shows that the mean relative error at test stations was generally about 20%. This implies that TEC predictions could recover about 80% of the observation TEC. It can be seen in Table 2 that during the two severe geomagnetic storm periods, the model relative error was degraded with respect to the 3-day mean value by a few percent to 10%. As stated previously, the relative error was only a ‘‘relative indicator’’ of assessing the ionospheric model prediction accuracy because of its dependence on the absolute value of the TEC observation itself. During periods I and II, the increased concentration of electrons in the ionosphere resulted in larger TEC observation values, as shown in Figure 4. This led to a smaller relative error, which did not increase as significantly as the VTEC RMS errors. At the IDTD station, the mean relative error for period II was

9 of 12


RS1007

RS1007

Figure 6. Slant TEC prediction relative error at test stations. actually lower than the 3-day mean value by 11.41%. Referring to Table 1, the IDTD station suffered a degradation of 1.937 TECU in the absolute VTEC prediction accuracy for period II, translating to a 37.15% accuracy loss compared to the 3-day mean accuracy. The results at IDTD further showed that the VTEC RMS error was more objective than relative

error measure in reflecting variations in ionospheric prediction accuracy.

5. Conclusions [16] This paper investigates the short-term (5-min) prediction performance of an ionosphere tomographic

Table 1. Statistics of Mean Vertical TEC Prediction RMS Errorsa Period I

AMC2 ERLA HRN1 IDTD PCK1

Period II

29 – 31 October 2003

RMS Error, TECU

Degradation, TECU

Degradation Percentage

RMS Error, TECU

Degradation, TECU

Degradation Percentage

RMS Error, TECU

7.32 7.46 7.48 7.57 7.70

2.60 2.48 1.39 2.36 1.15

55 50 23 45 18

7.94 8.06 8.47 7.15 8.65

3.22 3.08 2.38 1.94 2.10

68 62 39 37 32

4.72 4.98 6.09 5.21 6.55

a

Units are in TECU (1 TECU = 1016 el m – 2).

10 of 12


RS1007

RS1007

Table 2. Statistics of Mean Slant TEC Prediction Relative Errors

AMC2 ERLA HRN1 IDTD PCK1 Mean of five stations

Relative Error for Period I, %

Relative Error for Period II, %

Relative Error for 29 – 31 October 2003, %

34 10 32 26 17 24

22 10 31 9 16 17

26 8 28 20 16 20

model under high levels of ionospheric activity using GPS data recorded during the October 2003 geomagnetic storm. The mathematical formulation for tomographic modeling is developed, and two measures of prediction accuracy are computed. Three days of GPS data are analyzed, as collected from a wide-area GPS network including 55 stations in North America. Accuracy estimates are derived at five truth test sites, and results indicate that a mean prediction accuracy of 5.0 6.5 TECU can be achieved in the vertical direction. Results also show that the model’s prediction accuracies are degraded by 2.0 3.0 TECU for vertical TEC predictions during periods with extremely high ionospheric activities (Kp = 9). This is equivalent to a 38% 48% accuracy loss in comparison to the 3-day mean accuracy. Another accuracy assessment indicator, relative error, shows that the TEC predictions could statistically recover approximately 80% of observed TEC. The relative error statistics also show that TEC predictions are degraded by several percent to 10% from the average performance during the two severe storm periods with Kp = 9.

References Bisnath, S. B. (2000), Efficient, automated cycle-slip correction of dual-frequency kinematic GPS data, paper presented at ION GPS 2000, the 13th International Technical Meeting of the Satellite Division of the Institute of Navigation, Salt Lake City, Utah, 19 – 22 Sept. Blewitt, G. (1990), An automatic editing algorithm for GPS data, Geophys. Res. Lett., 17(3), 199 – 202. Brown, R. G. (1983), Introduction to Random Signal Analysis and Kalman Filtering, 347 pp., John Wiley, Hoboken, N. J. Coster, A. J., J. C. Foster, and P. J. Erickson (2003), Monitoring the ionosphere with GPS, GPS World, 14(5), 42 – 49. Daniell, R. E., Jr., L. D. Brown, D. N. Anderson, M. W. Fox, P. H. Doherty, D. T. Decker, J. J. Sojka, and R. W. Schunk (1995), Parameterized ionospheric model: A global ionospheric parameterization based on first principles models, Radio Sci., 30(5), 1499 – 1510.

Department of Commerce (2004), NOAA Service assessment: Intense space weather storms, October 19 – November 07, 2003, Silver Spring, Md. Doherty, P., A. J. Coster, and W. Murtagh (2004), Space weather effects of October – November 2003, GPS Solutions, 8, 267 – 271. Flores, A., G. Ruffini, and A. Rius (2000), 4D tropospheric tomography using GPS slant wet delays, Ann. Geophys., 18, 223 – 234. Hajj, G. A., and L. J. Romans (1998), Ionospheric electron density profiles obtained with the Global Positioning System: Results from the GPS/MET experiment, Radio Sci., 33(1), 175 – 190. Hansen, A. J., T. Walter, and P. Enge (1997), Ionospheric correction using tomography, paper presented at ION GPS 1997, the 10th International Technical Meeting of the Satellite Division of the Institute of Navigation, Kansas City, Mo., 16 – 19 Sept. Howe, B. M., K. Runciman, and J. A. Secan (1998), Tomography of the ionosphere: Four-dimensional simulations, Radio Sci., 33(1), 109 – 128. Klobuchar, J. A. (1986), Design and characteristics of the GPS ionospheric time delay algorithm for single-frequency users, paper presented at PLANS-86 Conference, Inst. of Electr. and Electron. Eng., Las Vegas, Nev., 4 – 7 Nov. Klobuchar, J. A. (1987), Ionospheric time-delay algorithm for single-frequency GPS users, IEEE Trans. Aerosp. Electron. Syst., 23(3), 325 – 331. Komjathy, A., L. Sparks, T. Mannucci, and A. Coster (2004), The ionospheric impact of the October 2003 storm event on WAAS, paper presented at ION GNSS 2004, the 17th International Technical Meeting of the Satellite Division of the Institute of Navigation, Long Beach, Calif., 21 – 24 Sept. Liu, Z. Z. (2004), Ionosphere tomographic modeling and applications using Global Positioning System (GPS) measurements, Ph.D. thesis, Dep. of Geomatics Eng., Univ. of Calgary, Calgary, Alberta, Canada. Liu, Z. Z., and Y. Gao (2003), Ionospheric TEC predictions over a local area GPS reference network, GPS Solutions, 8, 23 – 29. Liu, Z. Z., and Y. Gao (2004), Development and evaluation of a new 3D ionospheric modeling method, Navigation, 51(4), 311 – 329.

11 of 12

RS1007


Raymund, T. D., J. R. Austen, S. J. Franke, C. H. Liu, J. A. Klobuchar, and J. Stalker (1990), Application of computerized tomography to the investigation of ionospheric structures, Radio Sci., 25(5), 771 – 789. Skone, S. (1998), Wide area ionosphere grid modeling in the auroral region, Ph.D. thesis, Dep. of Geomatics Eng., Univ. of Calgary, Calgary, Alberta, Canada. Skone, S. H. (2001), The impact of magnetic storms on GPS receiver performance, J. Geod., 75(9 – 10), 457 – 468. Woods, T. N., F. G. Eparvier, J. Fontenla, J. Harder, G. Kopp, W. E. McClintock, G. Rottman, B. Smiley, and M. Snow

RS1007

(2004), Solar irradiance variability during the October 2003 solar storm period, Geophys. Res. Lett., 31, L10802, doi:10.1029/2004GL019571.

Y. Gao and S. Skone, Department of Geomatics Engineering, University of Calgary, 2500 University Drive N.W., Calgary, AB, Canada T2N 1N4. Z. Liu, Leica Geosystems, Ltd., Suite 112, 5438 11th Street N.E., Calgary, AB, Canada T2E 7E9. ([email protected])

12 of 12