Performance Evaluation of Ionosphere Modeling Using Spherical Harmonics in the Korean Peninsula

Abstract

The signal broadcast from a GPS satellite experiences code delay and carrier phase advance while passing through the ionosphere, which causes a signal error. Many ionosphere models have been studied to correct this ionospheric delay error. In this paper, the ionosphere modeling for the Korean Peninsula was carried out using a spherical harmonics based model. In contrast to the previous studies, we considered a real-time ionospheric delay correction model using fewer number of basis functions. The modeling performance was evaluated by comparing with a grid model. Total number of basis functions was set to be identical to the number of grid points in the grid model. The performance test was conducted using the GPS measurements collected from 5 reference stations during 24 hours. In the test result, the modeling residual error was smaller than that of the existing grid model. However, when the number of measurements was small and the measurements were not evenly distributed, the overall trend was found to be problematic. For improving this problem, we implemented the modeling with additional virtual measurements.

1. INTRODUCTION

The global positioning system (GPS) refers to a system which can precisely trace a user position using the signal from the satellite, and is currently utilized in various fields. For the signal from the satellite, a delay occurs before reaching the user while passing through the ionosphere and troposphere, which causes an error in the user position. Therefore, diverse models have been developed to correct this error. The ionospheric delay error changes significantly by the solar activity, and has a large variation depending on location. For the last 20 years, many studies have been performed to precisely model the ionosphere which has the above-mentioned characteristics, and the studies are also currently in progress. Depending on purpose, the ionosphere model can be divided into the real-time model for providing the correction information to user in real time, and the post-processing model for various research purposes. The representative real-time models are the Klobuchar model which is currently broadcast from GPS satellites by being included in the navigation data, and the grid model which is used for the Wide Area Differential GPS (WADGPS) such as the Wide Area Augmentation System of the United States. The post-processing models are the Global Ionospheric Map model which provides the latitude and longitude of the entire globe by dividing into a grid of 2.5° × 5°, and the International Reference Ionospheremodel (Choi et al. 2010).

In this paper, the ionosphere modeling for the Korean Peninsula was carried out using spherical harmonics. A number of studies have already been performed regarding the ionosphere modeling using spherical harmonics (Choi et al. 2010, Moon 2004, Venkata Ratnam et al. 2009). However, most previous studies have used spherical harmonics for the post-processing model, and the performance has not been evaluated when spherical harmonics are used for the real-time ionospheric correction model. In this paper, the spherical harmonics ionosphere modeling was performed which appropriately adjusted the number of basis functions considering the real-time ionospheric correction, and assuming that the model is used for WADGPS, the performance was tested by comparing with the grid model which is currently in use. The actual data collection and Differential Code Bias (DCB) estimation procedure were conducted in a post-processing way, and the modeling was performed using the obtained ionospheric delay measurement.

2. IONOSPHERIC DELAY MEASUREMENT

The signal transmitted from a satellite passes through the atmosphere before reaching the user. In the upper atmosphere, gas molecules are ionized by the effect of ultraviolet rays and exist as the plasma state which has a high density of free electrons. This region is called ionosphere, and the propagation speed of a signal becomes slower when it passes through the ionosphere, which causes an error. The ionospheric delay error is proportional to the distance that the satellite signal travels before reaching the user, and thus the error gets larger as the satellite elevation angle decreases.

Fig. 1. Single layer ionosphere model (Moon 2004).

The ionosphere is distributed over a wide altitude range. However, when estimating the ionospheric delay value, it is generally assumed that the electron density is concentrated at specific altitude as shown in Fig. 1, and the method which multiplies the slope coefficient is used. The slope coefficient is a function of satellite elevation angle, and as shown in Eq. (1), the ionospheric delay value \((I_S)\) is expressed as the product of the vertical ionospheric delay value \((I_V)\) and the slope coefficient (Kim 2007).

\(I_S = I_V \cdot Q = I_V \sqrt{1- \left ( {\cos E \over 1+ {H \over R} } \right )^2}\)                                                                          (1)

\(E\): Satellite elevation angle
\(H\): Single layer ionosphere altitude (350 km)
\(R\): Earth’s radius
\(Q\): Slope coefficient

One of the important characteristics of ionospheric delay is that the degree of delay and advance for the code and carrier of a signal caused by the ionosphere changes depending on the frequency of the signal. With this characteristic, a dual frequency user can obtain the ionospheric delay value using the L1 frequency measurement and L2 frequency measurement. When the ionospheric delay is calculated using the difference of L1 and L2 frequency pseudoranges, the measurement noise of pseudorange is involved, and the relevant error occurs in the estimated ionospheric delay value. Therefore, this needs to be eliminated for obtaining a more precise ionospheric delay value. The Hatch filter can remove the noise of pseudorange measurement using the carrier phase measurement which has much smaller measurement noise than the pseudorange. The weighted Hatch filter combines the Hatch filter and the weighted least square method. For the weighted Hatch filter, the measurement of current epoch is estimated by assigning weight, which considers each noise variance value, to the estimated measurements of up to the previous epoch and the current measurement. In this regard, the standard deviation of measurement noise was modeled as an exponential function of elevation angle for each reference station, and the filter equations are shown in Eqs. (2-4) (Han et al. 2012).

\( \widehat{d_{\rho_k}} = { \hat{\sigma}_{d\rho,k}^2 \over \hat{\sigma}_{d\rho,k-1}^2} [ \hat{d\rho}_{k-1} - \Delta(d\phi_k) ] + { \hat{\sigma}_{d\rho,k}^2 \over \sigma_{d\rho,k}^2} d\rho_k \)                                                                 (2)

\( {1 \over \hat{\sigma}_{d\rho,k}^2} = {1 \over \hat{\sigma}_{d\rho,k-1}^2} + {1 \over \sigma_{d\rho,k}^2} \)                                                                                       (3)

\(\hat{i}_k = { \widehat{d\rho_k} \over \gamma -1}\)                                                                                             (4)

where

\( d_{ρk}\): L1 and L2 pseudorange difference measurement \((𝜌_{2,k}−𝜌_{1,k)}\) at \(k\)-th epoch
\( \hat{d_{ρk}}\): L1 and L2 pseudorange difference filter estimate at k-th epoch
\( \sigma_{d\rho,k}\): Noise standard deviation of \(d\rho_k\) at \(k\)-th epoch
\(\hat{\sigma}_{d\rho,k}\): Noise standard deviation of \(d\rho_k\) at \(k\)-th epoch
\(d \phi_k\): L1 and L2 carrier difference measurement \((ϕ_{2,k}−ϕ_{1,k)}\) at \(k\)-th epoch
\(\Delta(d\phi_k)\): Difference between the \(d \phi\) value of \(k\)-th epoch and the \(d \phi\) value of \((k-1)\)-th epoch
\(i_k\): Estimated ionospheric delay value of \(k\)-th epoch
\(\gamma\): Ratio of the square of L1 and L2 frequency \((\approx 1.65)\)

3. IONOSPHERE MODELING

3.1 Grid Model

The grid model estimates and provides the predefined vertical ionospheric delay value of a grid point, and the ionospheric delay is corrected by calculating the ionospheric delay value at a specific ionospheric pierce point from the ionospheric delay value of the grid point. For the standard grid model of WADGPS, the grid interval of 5° in the latitude and longitude directions is generally used. Fig. 2 shows the conceptual diagram for estimating the vertical ionospheric delay value of a grid point from the ionospheric delay measurements, and the ionospheric delay value of the grid point is estimated using the measurements that are within a certain radius of each grid point. The weighted interpolation is used as the estimation method, and the calculation equation is shown in Eq. (5) (Kim 2007). The ratio of ionospheric delay measurement for the Klobuchar model value is interpolated rather than simply interpolating the ionospheric delay measurement, and this is to secure the trend for the local time.

\( I_{Grid,V} = I_{Klob}^{Grid} { \sum\limits_{k=1}^K [ \left ( I_{meas,V}^k \over I_{Klob,V} \right ) \over \sum\limits_{k=1}^K w_k } \)                                                                               (5)

where

\(I_{Grid,V}\): Vertical ionospheric delay estimated at a grid point
\(I_{Klob}^{Grid}\): Vertical ionospheric delay estimated at a grid point using Klobuchar model
\(I_{Klob,V}^k\): Vertical ionospheric delay estimated at the ionospheric pierce point of a measurement using Klobuchar model
\(I_{meas,V}^k\): Measured vertical ionospheric delay
\(\omega_k\): Weight for each measurement

Fig. 2. Conceptual figure of grid point and pierce point.

The weight for each measurement is a value that can vary depending on how it is defined by the user. In this paper, the weight was made to be inversely proportional to the distance and measurement variance. The ionospheric delay value at the ionospheric pierce point, which is to be obtained from the ionospheric delay value of the grid point, is calculated in a similar way to the above procedure.

3.2 Spherical Harmonics Model

The spherical harmonics refer to the basis functions which are orthogonal one another on a sphere as shown in Fig. 3. A certain data on a sphere can be expressed as a linear combination of these basis functions, and when the vertical ionospheric delay value is expressed as spherical harmonics, it is as shown in Eq. (6) (Kim 2007).

\(I_V(\phi, \lambda)=\sum\limits_{n=0}^N \sum\limits_{m=0}^N P_{nm} (\sin \phi) (C_{nm} \cos (m\lambda) +S_{nm} \sin (m\lambda))\)                                                     (6)

where

\(I_V\): Vertical ionospheric delay
\(\phi\): Magnetic latitude of an ionospheric pierce point
\(\lambda\): Local time of an ionospheric pierce point which normalized 24 hours as \(2\pi\)
\(N, M\): Maximum degree and order
\(C_{nm},S_{nm}\): Spherical harmonics coefficient
\(P_{nm}\) : associated Legendre function

Fig. 3. Spherical Harmonic functions.

In Eq. (6), the vertical ionospheric delay value is obtained from the L1 and L2 frequency measurements, and the magnetic latitude and local time of an ionospheric pierce point are calculated from the latitude and longitude for the ionospheric pierce point of a measurement. n and m refer to degree and order, respectively. The degree determines resolution in the latitude direction, and the order determines resolution in the longitude direction. In other words, as the maximum degree and order values become larger, more detailed and precise expressions are possible in the latitude and longitude directions, respectively. In theory, as the maximum degree and order increase, more precise modeling is enabled, and if it is used as a post-processing model for the entire globe, the values are typically set to more than 10. However, if it is used as a model for providing the real-time correction, the maximum degree and order values need to be adjusted appropriately considering computation and the delivery of correction. The spherical harmonics coefficient \(C_{nm}, S_{nm}\) is the unknown which is to be estimated. If the spherical harmonics expression for the measurement is rearranged by letting this coefficient be x, it is as shown in Eq. (7).

\(z= \left [ \begin{matrix} I_{V,1} \\ I_{V,2} \\ \vdots \\ I_{V,k} \end{matrix} \right ] = H_x\)                                                                                 (7)

where

\(H_{i,j}=\left[\begin{matrix}P_{ij}(\sin{\phi_1})\cdot \cos{(}j\cdot\lambda_1)&P_{ij}(\sin{\phi_1})\cdot \sin{(}j\cdot\lambda_1)\\P_{ij}(\sin{\phi_2})\cdot \cos{(}j\cdot\lambda_2)&P_{ij}(\sin{\phi_2})\cdot \sin{(}j\cdot\lambda_2)\\\vdots&\vdots\\P_{ij}(\sin{\phi_k})\cdot \cos{(}j\cdot\lambda_k)&P_{ij}(\sin{\phi_k})\cdot \sin{(}j\cdot\lambda_k)\\\end{matrix}\right]\)

\(H=\left[\begin{matrix}H_{0,0}&H_{0,1}&\cdots&H_{N,M}\\\end{matrix}\right]\)

\(x=\left[\begin{matrix}C_{00}&S_{00}&\cdots&\begin{matrix}C_{NM}&S_{NM}\\\end{matrix}\\\end{matrix}\right]\)

From Eq. (7), the coefficient is estimated as shown in Eq. (8) using the least square method.

\( x = (H^TH)^{-1}H^Tz \)                                                                                   (8)

When the coefficient is obtained, the vertical ionospheric delay value at a specific ionospheric pierce point can be calculated from Eq. (6).

4. MODELING PERFORMANCE TEST

4.1 Data Collection

The performance test was carried out using the measurements obtained from the NDGPS reference stations of the Ministry of Land, Transport and Maritime Affairs. The Daejeon, Palmido, Eocheongdo, Sokcho, and Dokdo reference stations were used, and a day (July 2, 2012) was selected which has no anomalies such as ionospheric storm. The ionospheric delay measurement was obtained by filtering the data which was stored every second in the BINary EXchange format and the data for a 24-hour period was used.

When the data from a number of reference stations are used, the DCB value should be removed which is induced by the electrical path difference between the L1 frequency signal processing circuit and the L2 signal processing circuit. Every satellite DCB value can be removed in the data processing procedure of each reference station using the Time of Group Delay value of the navigation data provided by the GPS satellite. The DCB value for each reference station is calculated by collecting the data from every reference station, and can be estimated along with the spherical harmonics coefficient by adding the reference station DCB term to the spherical harmonics expression of the ionospheric delay measurement (Kim 2007). The test was conducted in a post-processing way. As the reference station DCB value does not change greatly for a day, the DCB value was estimated by processing the data for the entire test period at a time, and the estimated DCB value was removed from every ionospheric delay measurement prior to using the measurement. Table 1 lists the estimated reference station DCB values.

Table 1. DCB of reference station.

Reference station	DCB (m)
Daejeon Eocheongdo Sokcho Dokdo Palmido	-4.6247 -6.2198 -4.8157 -5.9887 -5.9101

 

4.2 Performance Comparison Between Spherical Harmonics Model and Grid Model

For the grid model, the modeling was performed using total 16 grid points as shown in Fig. 4. For the spherical harmonics model, the maximum degree and order values were selected so that the number of spherical harmonics basis functions also becomes 16. In this case, both the maximum degree and order values were 3.

Fig. 4. Position of Korea NDGPS reference stations (blue) and Grid Point (red).

To compare the modeling performance of the two models, the modeling residual error was calculated at each time, and Fig. 5 shows the result. The modeling residual error indicates the difference between the measurement value used for modeling and the ionospheric delay value estimated at the ionospheric pierce point of the measurement for each model. The top figure is a graph that shows the maximum value among the residual errors at each time, and the bottom figure is a graph that shows the rms value of the residual errors.

Fig. 5. Time history of residual error.

As shown in Fig. 5, the modeling residual error was similar between the spherical harmonics model and the grid model, but the spherical harmonics model showed slightly smaller values. On average, the maximum residual error of the spherical harmonics model was lower by 0.24 m, and the rms value of the spherical harmonics model was lower by 0.09 m.

To compare the overall shape of ionosphere, the result for each model was presented using a figure. Fig. 6 shows the vertical ionospheric delay values of the two models at the GPS week seconds 93,630. The left figure is the result of the spherical harmonics model, and the right figure is the result of the grid model. The pink dot represents the ionospheric pierce point of the measurement. Fig. 7 shows the result which processed the data of IONosphere Exchange (IONEX) format at the corresponding time. The left figure is the result for the entire globe, and the right figure is the result for the Korean Peninsula. The IONEX data is provided having a grid of 2.5°×5° (latitude and longitude, respectively) for the entire ionosphere of the globe. As the measurement used for the modeling is different from the IONEX data, the modeling result cannot be directly compared with the result of IONEX. However, it was presented to be used as the criterion for the approximate value and shape.

Fig. 6. The modeling result of vertical Ionospheric delay at GPST=93630 (Spherical harmonics, grid model).

Fig. 7. The IONEX data at GPST=93630.

As for the overall ionospheric delay value of the Korean Peninsula shown in Fig. 7, the value decreases with increasing latitude, and the value increases with increasing longitude. Both spherical harmonics model and grid model were noticeably different from the IONEX data in the low latitudes near 135° longitude. It is not thought that this is solely due to the error from modeling because the ionospheric delay measurement, which was measured at the corresponding region, was different from the IONEX model. However, the spherical harmonics model showed a considerably different trend from the IONEX data on every edge and in the region without measurement. Especially, in the vicinity of 45° N, 120° E, the value abruptly decreased, and in the vicinity of 30° N, 120° E, the value should decrease with decreasing longitude, but the result showed that the value increases. As the grid model estimates the result by interpolating the data, the trend of data is maintained even in the region without measurement. In contrast, the spherical harmonics model estimates the coefficient without considering the overall trend, and tries to similarly match the value only at the measurement location, which leads to the above-mentioned results. In other words, if the number of measurements is small and the measurements are not evenly distributed, the spherical harmonics model is not able to estimate the value reflecting the overall trend in the region without measurement, which makes it hard to be used as a correction model.

4.3 Performance Evaluation of Spherical Harmonics Model Using Virtual Measurements

If the previous result is considered the other way around, as the spherical harmonics model had a smaller modeling residual error for the measurement, the modeling result reflecting the overall trend of the measurement can be obtained when measurements are evenly distributed. Therefore, in this paper, the modeling was performed by adding virtual measurements on the edge of the region. For the virtual measurements, the 12 grid points on the edge, estimated from the grid model, were used, and in Fig. 8, the orange dot represents the point. When estimating the coefficient, to prevent the deterioration of modeling performance for the actual measurement, the weight of the added virtual measurement was set to 1/10 of the weight of the actual measurement. To differently assign the weight in this manner, Eq. (8) for estimating the coefficient was converted to the weighted least square method as shown in Eq. (9).

\(x=(H^TWH)^{-1}H^TW_z\)                                                                           (9)

\(W = \left [ \begin{matrix} 1 & \cdots & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & 0 & 0 & \cdots & 0 \\ 0 & 0 & 1 & 0 & \cdots & 0 \\ 0&0&0&1/10&\cdots&0 \\ \vdots &\vdots &\vdots &\vdots &\ddots &\vdots \\ 0&0&0&0&\cdots&1/10 \end{matrix} \right ]\)

where, \(W\): Weight for the measurement.

Fig. 8. The modeling result of vertical Ionospheric delay at GPST=93630 (Spherical harmonics, grid model).

The modeling residual error was also calculated when the virtual measurements were used, and Fig. 9 shows the overall diagram for the modeling result. For the spherical harmonics model, when the virtual measurements were added, the modeling residual error for the measurement was rather increased compared to when the virtual measurements were not added, but the residual error was still slightly smaller than the grid model. On average, the maximum residual error was smaller by 0.21 m, and the rms value was smaller by 0.05 m. In Fig. 8, the overall trend of ionosphere, which was previously problematic, is found to be to some degree improved. The region that showed a considerably different trend from the actual trend of ionosphere near 120° longitude when the virtual measurements were not added was improved, but in the vicinity of 40°~45° latitude, the expected trend, where the value decreases with increasing latitude, was still not observed. Based on the above result, it is thought that the overall trend of the value can be improved without significantly increasing the modeling residual error for the existing actual measurement when the virtual measurement, to which relatively small weight is assigned, is utilized.

Fig. 9. Time history of residual error.

5. CONCLUSION

In this paper, the ionosphere modeling for the Korean Peninsula was carried out using spherical harmonics. The modeling performance was evaluated by comparing with the grid model. The number of basis functions was set to be identical to the total number of grid points in the grid model considering real-time modeling. The modeling was performed for each model using the measurements obtained from 5 NDGPS reference stations. To evaluate the performance, each modeling residual error was calculated and the overall shape was presented using a figure. In the analysis of the size of residual error, the spherical harmonics model had a smaller residual error than the grid model. However, for the overall shape which was compared with the IONEX data, the spherical harmonics model was not able to properly express the trend depending on the changes of latitude and longitude in the region without measurement, which makes it hard to be used as a real-time correction model as it stands. To improve this, the values of grid points on the edge of the region, which were estimated from the grid model, were used as the virtual measurements. By applying the virtual measurements to the spherical harmonics model, it was found that the overall modeling shape can be improved without significantly increasing the modeling residual error for the existing measurement. However, as the suggested result is based on the data for a 24-hour period, the test and analysis should be carried out for a longer period of time. Also, the generation method and location of virtual measurements need to be studied to secure the overall trend properly.