Global GPS Ionospheric Modelling Using Spherical Harmonic Expansion Approach
 Author: Choi ByungKyu, Lee Woo Kyung, Cho SungKi, Park JongUk, Park PilHo
 Organization: Choi ByungKyu; Lee Woo Kyung; Cho SungKi; Park JongUk; Park PilHo
 Publish: Journal of Astronomy and Space Sciences Volume 27, Issue4, p359~366, 15 Dec 2010

ABSTRACT
In this study, we developed a global ionosphere model based on measurements from a worldwide network of global positioning system (GPS). The total number of the international GPS reference stations for development of ionospheric model is about 100 and the spherical harmonic expansion approach as a mathematical method was used. In order to produce the ionospheric total electron content (TEC) based on grid form, we defined spatial resolution of 2.0 degree and 5.0 degree in latitude and longitude, respectively. Twodimensional TEC maps were constructed within the interval of one hour, and have a high temporal resolution compared to global ionosphere maps which are produced by several analysis centers. As a result, we could detect the sudden increase of TEC by processing GPS observables on 29 October,2003 when the massive solar flare took place.

KEYWORD
global positioning system , global ionosphere model , spherical harmonic expansion , total electron content

1. INTRODUCTION
The modeling of the ionosphere using global positioning system (GPS) has been intensively studied in last two decades and the research is still actively performed. Since it changes rapidly in terms of time and space, modeling of the ionosphere is not easy. However, 32 GPS satellites and the geometrical arrangement of the GPS reference stations distributed on the ground allow the signal analysis of various paths. Since the ionosphere is the largest error source on the path through which the GPS signal is transmitted, it reduces the accuracy of the GPS user’s position and the navigation information. While the previous ionosphere models were developed as the means to reduce the delay error of GPS signal and improve the accuracy of the GPS user’s position, current models are focused on the detailed study and analysis of the scientific features of the ionosphere through the combination with the information received by various observation instruments such as ionosonde and incoherent scatter radar.
The GPS ionosphere models can be divided into the realtime models that are used in positioning and navigation and the posttime models that can be used for various scientific purposes. The former category includes the ionosphere grid models that are applied to the global navigation augmentation system such as Klobuchar model and wide area argumentation system (WAAS) /European geostationary navigation overlay service (EGNOS)and the latter category is represented by global ionosphere maps (GIM) and global assimilation ionospheric measurements (GAIM) (HofmannWellenhof et al. 2008). The GPS ionosphere models are also divided into the local models confined in a specific location and the global models. The local models have high temporal and spatial resolution using the observation information received from the GPS reference stations that are densely distributed (Mannucci et al. 1998). On the contrary,the global models employ the measurements from the international GNSS service (IGS) GPS reference stations that are internationally authenticated. Currently, there are about 420 international GPS reference stations registered to the IGS and about 100 of them are used in the global ionosphere modeling. This number depends on the ionosphere modeling strategies of the international ionosphere analysis centers (ACs) (Mannucci et al. 1998, HernandezPajares et al. 1999). In general, the ionosphere ACs select the GPS reference stations appropriate to the ionosphere modeling and the foremost consideration of the selection criteria is the geographical positions of the GPS reference stations distributed on the ground. The institutes that serve as ionosphere ACs include Jet Propulsion Laboratory (JPL) in USA, Center for Orbit Determination in Europe (CODE) in Switzerland, Natural Resource Canada (NRCan) in Canada, European Space Operation Center (ESOC) in Germany and Polytechnical University of Catalunya (UPC) in Spain. These institutes have developed independent global GPS ionosphere models and provide vertical total electron contents (TEC) values in the form of ionosphere exchange (IONEX) (Schaer 1999). The provided TEC values have the spatial resolution of 2.5° and 5° in latitude and longitude, respectively, and the temporal resolution of 2 hours. The global ionosphere model is the technology held by only a few advanced institutes and it involves processing of huge amount of data, long processing hours and constitution of difficult mathematical model. In this study, a global ionosphere model was developed applying the measurements from 100 international GPS reference stations distributed on the ground and the spherical harmonic expansion technique. We present the model that has relatively higher temporal resolution than that of the conventional models to quickly detect the effect on the geomagnetic storm on the global ionosphere.
2. ESTIMATION OF IONOSPHERE TEC
The ionosphere as an electrolyte medium, it causes GPS signal delay or carrier phase advance when electromagnetic wave like navigation signal passes through it. In addition, the signal delay is converted to the distance error as shown in Eq. (1):
where
f denotes the frequency of the electromagnetic wave, the carrier phase has positive (+) value and the pseudorange measurement has negative () value.Using the pseudorange measurements of GPS L1 and L2 simultaneously, the TEC in Eq. (1) is easily calculated by Eq. (2) (Blewitt 1990):
where
P_{i(i=1,2)} denotes the pseudorange measurements of GPS L1 and L2, respectively, andf_{i(i=1,2)} denotes the frequencies of L1 and L2. Each of the pseudorange measurements of Eq. (2) includes hardware bias. If the ‘GeometryFree Linear Combination’ is applied as shown in the Eq. (2) above, the hardware bias has the form of difference, which is referred as differential code bias (DCB) (Lanyi & Roth 1988, Coco et al. 1991). The DCB values are separated into the GPS satellite and the receiver and the magnitudes of the values are shown in Figs. 1 and 2. It is not easy to estimate the DCB values of the satellite and the receiver accurately from the GPS observations. Fortunately, the ionosphere ACs provides the DCB values of the GPS reference stations registered to the IGS and the GPS satellites. The DCB values in Figs. 1and 2 are the values provided by the IGS. The DCB values of the GPS satellites range between 3 ns and +4 ns and are very stable. The DCB values of the GPS reference stations were between the minimum of about 18 ns and the maximum 20 ns. In addition, it was found that the individual GPS satellites and the receivers had different DCB values. Comparing the satellite DCB and the receiver DCB, the receiver DCB was larger than that of the satellite in general, which means that the receiver DCB values have more effect on the TEC estimation than the satellite DCB values (Choi et al. 2010).
If the DCB values of the GPS satellites and receivers are determined, the above Eq. (2) can be rewritten as Eq. (3) (Blewitt 1990):
where
b_{r} andb^{s} denote the DCB values of receiver and satellite, respectively.3. MATHEMATICAL MODELING
The ionosphere is a region between 50 km and 1,000 km from the ground level where free electrons are densely distributed. Generally, it is assumed that free electrons in the ionosphere are concentrated at a certain altitude like a thin layer to for the 2dimensional modeling of the ionosphere TEC change. Many researchers choose the altitude between 250 km and 450 km (Tseng et al. 2010). In this study, we chose 350 km which is the altitude that is generally used in GIM model.
After determining the altitude of the electron concentration and the satellite, the mapping function can be simply expressed as Eq. (4):
where
z denotes the zenith angle of the reference station andz' denotes the zenith angle of the single layer. In addition, R is the radius of the earth andH is the height of the single layer from the ground (assumed to be 350 km).It is not an easy task to estimate largescale TEC such as global TEC from the GPS observations. Since the ocean takes approximately 70% of the earth surface, the installed GPS reference stations are concentrated on the land. Because of this, the observation data received from the land and the ocean are seriously biased and even the international GPS reference stations are concentrated particularly in North America and Europe regions. Such like distribution of the GPS reference station is directly related to the temporal and spatial resolutions of observation and mathematical model and causes reduced accuracy of the TEC estimation.
Most of the ionosphere ACs including JPL and CODE apply the spherical harmonic expansion technique to develop global ionosphere models, while Polytechnical University of Catalunya (UPC) uses the Kriging method based on interpolation (HernandezPajares et al. 1999). The advantage of the spherical harmonic expansion technique is that data of large capacity can be approximated by a few numbers of coefficients, while the disadvantage of it is that the information that varies drastically in a very short period of time is eliminated. Although the Kriging method which is the interpolating method where weight is applied depending on the distance is too sensitive to the measurements, it has the advantage that it can be conveniently applied to the detection of momentous change in a short period of time. The Kriging method is more appropriate to local models.
In this study, the spherical harmonic expansion technique which is generally used and convenient for the processing of large capacity data was applied to calculate the global TEC. The applied spherical harmonic expansion technique is shown in Eq. (5) (Schaer 1999):
where
n_{max} denotes the maximum degree of spherical harmonic expansion, P^{~}_{nm}=Λ(n.m )？ P_{nm} is the normalized associated Legendre function with the degree,n , and the order,m, anda_{nm} andb_{nm} are the TEC unknown coefficients of the spherical harmonics. In addition,β means the geocentric latitude and s is the sunfixed longitude of the point where the ionosphere penetrates. Eq. (6) can be used to calculated s:where
l_{o} denotes the geographic longitude. The normal function Λ is defined as in Eq. (7):where δ denotes the Kronecker delta.
In summary, to calculate the gridbased ionosphere TEC corresponding to the spatial resolution suggested in this study, the vertical TEC is firstly calculated using Eqs. (3) and (4) and then the spherical harmonic technique of Eq. (5) is applied.
4. PROCESSING OF THE GPS OBSERVATION DATA AND THE RESULTS
To develop global GPS ionosphere model, the observation data should be firstly obtained from the international GPS reference points that are uniformly distributed on the earth ground. As mentioned in the introduction, there are about 420 reference stations that serve as the international GPS reference points. Although the observation data from all these GPS reference stations can be used in the model development, the number can be reduced to about 150 by considering the geographical distribution and excluding the GPS reference stations that are close to each other. If the data receive ratio and data loss are considered additionally, the observation data from 80120 GPS reference stations are available. In this study, the IONEX file was downloaded from the IGS data center to obtain the observation information, the navi
navigation data and the DCB values of the GPS satellites and receivers of 100 GPS reference stations, as shown in Fig. 3. Since the size of the GPS observation data is large, preprocessing was performed to save the observation and navigation data for each observation hour in the form of metadata for the convenience of further data processing. In addition, the DCB values were reestablished for the individual GPS reference stations. Following the preprocessing procedure, time synchronization of the observation information of the GPS reference stations was carried out. Considering the two codes information (C1 or P1 and P2) and two carrier phases information (L1 and L2) received by one reference station, the size of the array observed during one day is about 100,000 × 4. Additionally considering the number of reference stations which is 100, the entire size is 100 × 100,000 × 4 that requires processing of huge amount of observation data. Determining the position of the GPS satellite (X, Y and Z) using the time information and the navigation data after the time synchronization, the elevation angle and azimuth angle between the individual reference stations and the satellites were calculated. Subsequently, the ionosphere pierce point which is the point at which the GPS signal passes through the ionosphere was calculated using the radius of the earth (~6,370 km) as well as the altitude which was assumed to be the altitude at which the electrons of the ionosphere are concentrated. After the visual direction TEC was determined based on the observation data, the vertical TEC could be determined using the mapping function which is expressed as a function of the satellite altitude. Finally, the spherical harmonic expansion technique was applied to calculated the vertical TEC corresponding to the grid of which latitude and longitude are 2.5° and 5°, respectively. A total of 100 coefficients were obtained by expanding the order and degree to 10, respectively.
Fig. 4 shows the map of the international GPS reference stations used in the global ionosphere model. Since the total number of the ionosphere cannot be accurately estimated even by a mathematical model in the case where the GPS observation data is not available, the GPS reference stations were selected so that they could be distributed uniformly on the earth ground considering the geographical factors. Fig. 5 shows the point at which the GPS signal passes through the ionosphere and the points of penetration in one hour were indicated together in the interval of 15 minutes.
To verify the accuracy of the global GPS ionosphere model developed in this study, it was compared with the GIM provided by the international ionosphere ACs. Be
cause the GIM models provide the TEC values in 2hour interval, as described in the introduction, the GPS ionosphere model of Korea Astronomy and Space Science Institute (KASI) should also generate the ionosphere TEC values in 2hour interval for the comparison and verification. Fig. 6 shows the comparison between the GIM of the ionosphere ACs and the KASIGIM model. The calculated values of the KASIGIM model were very similar to those of JPL, IGS and CODE, but they were relatively larger than those of European Space Agency (ESA). Table 1 is the list of the root mean square of the difference between values of KASIGIM model and the AC model. The AC presents the error of GIM as to be approximately 2~8 TEC unit (TECU) and the KASIGIM model calculated the ionosphere TEC within this error range.
Fig. 7 is the GIM calculated by processing the data received by the 100 GPS reference station on October 29, 2003. The twodimensional ionosphere map was calculated in the one hour interval. As seen in Fig. 7, the ionosphere TEC drastically increased at Fig. 7u and later.
X17 flare and coronal mass ejection (CME) were observed on October 28, 2003 and X10 flare (maximum Xray flux of 10^{3} W/m^{2}) was generated at 20:49 UTC (coordinated universal time) on the next day, October 29, 2003 (Plunkett 2005). At that, the TEC calculated by GPS was as high as the instantaneous maximum of 200 TECU. Fig. 8 is the GIM produced at 21 (UT) on October 29, 2003. It shows that the ionosphere TEC was drastically increased in the ?southwest region of USA and the west region of Chile. This may be because the upper atmosphere of the earth was directly affected by the geomagnetic storm caused by the solar flare. Thus, it was demonstrated that the global GPS ionosphere developed in this study could be used to monitor the ionosphere variation by external effects such as the space environmental changes. Moreover, the developed model has the advantage in detecting momentous change of the ionosphere because it has better temporal resolution than that of the 2hour interval ionosphere TEC map provided by the international ionosphere ACs.
5. CONCLUSIONS
In this study, a global ionosphere model was developed using the observation information received from the 100 GPS reference stations that are uniformly distributed on the earth ground. The spherical harmonic expansion technique was used as the mathematical model for the global ionosphere model, expanding the degree and the order to 10. The ionosphere model developed in this study has the grid resolution of 2.5° and 5° of the latitude and longitude, respectively, and it generates a twodimensional ionosphere map in the interval of one hour. In contrast to the international ionosphere ACs that usually provide the TEC map in 2hour interval, our model generates the TEC map in the shorter temporal resolution which is the interval of one hour. There are two significant points of this model: the variation of the ionosphere can be quickly recognized by the increased temporal resolution; and the global ionosphere analysis is possible as KASI develops the technology independently although it is possessed by only a few international ACs. To observe the ionosphere variation characteristics using the developed model, we processed the GPS observation data of October 29, 2003 when the geomagnetic storm by the solar flare was detected. The result showed that the drastic increase of the ionosphere TEC was detected from the 21 hour (UT) on October 29, 2003. The detection of the instantaneous variation of the ionosphere has a great value in scientific research and other areas. Moreover, additional verification of the global ionosphere model developed in this study is required and we will perform the study to increase the spatial resolution considering the potential linkage to the local models in the future.

1. Blewitt G 1990 [GeoRL] Vol.7 P.199

2. Choi B. K, Chung J. K, Cho J. H 2010 [JASS] Vol.27 P.123

3. Coco D, Coker C, Dahlke S. R, Clynch. J. R 1991 [ITAES] Vol.27 P.931

4. HernandezPajares M, Juan J, Sanz J 1999 [JASTP] Vol.61 P.1237

6. Lanyi G. E, Roth T 1988 [RaSc] Vol.23 P.483

7. Mannucci A. J, Wilson B. D, Yuan D. N, Ho C. H, Lindqwister U. J, Runge T. F 1998 [RaSc] Vol.33 P.565

8. Plunkett S 2005 [NRL Review] P.91

9. Schaer S 1999

10. Tseng K. H, Shum C. K, Yi Y, Dai C, Lee H, Bilitza D, Komjathy A, Kuo C. Y, Ping J, Schmidt M 2010 [Mar Geodes] Vol.33 P.272

[Fig. 1.] A daily set of global positioning system satellite differential code bias values.

[Fig. 2.] A daily set of global positioning system receiver differential code bias values.

[Fig. 3.] Flow chart of data processing for global GPS ionosphere model.

[Fig. 4.] Global positioning system reference stations used for global ionosphere model (the marked circles on map mean the positions of them).

[Fig. 5.] Ionosphere pierce points calculated by using the measurements received from each global positioning system reference station.

[Fig. 6.] Comparison of the daily total electron content variations between KASI and analysis centers.

[Table 1.] Comparison of RMS values between KASI and ACs.

[Fig. 7.] Twodimensional ionospheric total electron content maps produced from GPS measurements on 29 October 2003. The first image (a) means the snap shot at 1:00 UT the last image (x) is at 24:00 UT.

[Fig. 8.] Global ionospheric total electron content map produced by KASI model at 21:00 UT on 29 October 2003.