Characterization of Ionospheric Delay and Forecasting Using GPS-TEC over Equatorial Region

This paper investigates the characteristic of the ionospheric delay and forecasting using modified statistical Holt-Winter method over Malaysia. This was carried out by using dual frequency GPS Ionospheric Scintillation and TEC Monitor (GISTM) receivers installed at the National University of Malaysia (UKM) (geographic coordinates 2.55°N, 101.46°E and geomagnetic coordinates 7.10°S, 174.05°E), and Langkawi national observatory (Langkawi) (geographical coordinates 6.19°N, 99.51°E and geomagnetic coordinates 3.39°S, 172.42°E), during 2011. In this work, the Holt-Winter forecasting method is modified to provide a better forecasting result over Malaysia. The diurnal, monthly, and seasonal variations of actual ionospheric delay have been analyzed and compared with forecasted ionospheric delay. Based on the results, the diurnal variability of ionospheric delay was maximal mostly between 12:00– 17:00 LT, and minimal nearly at 05:00 LT. The maximum monthly ionospheric delays occurred in October, which is attributed to the geomagnetic storm disturbance that caused increase in ionospheric delay. The lowest seasonal ionospheric delay occurred in the summer and the highest was during the equinox. The maximum ionospheric delay on a quiet day was around 41.9 TECU, while it reached up to 86.6 TECU during the geomagnetic storm. The Holt-Winter method exhibited the highest difference between the measurements and the forecasts in March, with an error of 6%. Hence, the modified Holt-Winter method forecast ionospheric delay effectively during quiet and disturbed periods over Malaysia..


Introduction
The Earth's ionosphere is the ionized region of the upper atmosphere consisting of different ionized layers (D, E, and F). The electron concentrations reach their highest values at F layer, which is largely affects the radio waves propagation. The Total Electron Content (TEC) is an important ionospheric parameter measuring the total number of electrons in a column of one square meter cross section along the path from satellite to receiver [Langley, 2000;Abdullah et al., 2009]. The ionospheric TEC is a notable source of error that disrupts the accuracy of the global positioning system (GPS) signals passing from the satellite to the ground by changing the velocity of the propagated signals causing ionospheric delay. The ionospheric delay changes with the consideration of factors including the user location, the elevation angle, the time of day, the seasons of the year, and the solar cycle. During high solar activity, the delay can generate vertical errors in the range of 5 and 15 m, and it can reach up to 100 m in particular conditions [Klobuchar, 1991]. In the low latitude and equatorial regions, the variations of the TEC are recognised to be high when compared with other regions such as mid-latitude and high-latitude. The ionospheric delay in equatorial regions is subject of day-to-day variability and is a challenging problem for ionospheric modelers [Kumar et al., 2016]. Therefore, investigation the ionospheric delay and identification of the precise and accurate forecasting model of transionospheric propagation errors are essential for precise measurements and further contributes valuable information to satellite and space probe navigation, space geodesy, radio astronomy and other applications.
A number of empirical models have been developed to effectively forecast ionospheric delay variability at various regions worldwide, such as the Bent model [Bent et al., 1975], Global Ionospheric Maps (GIMs) [Mannucci et al., 1998], and the International Reference Ionosphere (IRI) model [Bilitza, 2001], which is widely used in different scientific areas and continuously developed and updating. Generally, the ionospheric delay error correction for single frequency GPS receivers is done using the broadcast model, with the accuracy ranging from 50% to at most 60% [Klobuchar, 1996]. Further techniques have been developed for ionospheric prediction purposes such as the Autoregressive Moving Average (ARMA) method, and the autocovariance that is used in forecasting the TEC timeseries in various stations over Europe, demonstrating a precision forecast in the mid latitudes [Krankowski et al., 2005;Stanislawaka and Zbyszynski, 2002]. Furthermore, the autocovariance forecasting procedure was used to forecast ionospheric delay in disturbed conditions with acceptable precision [Stanislawska and Zbyszynsk, 2001;Stanislawska and Zbyszynsk, 2002]. A number of researchers have used self-consistent and neural network models for forecasting the TEC over the equatorial region [Williscroft and Pool, 1996;Oyeyemi et al., 2006;Habarulema et al., 2009, Klimenko et al., 2011. Later, Suwantragul and other researchers demonstrated the application of the Holt-Winter method in forecasting the ionospheric delay over Chiang Mai in Thailand using the GPS-TEC measurements [Suwantragul et al., 2003]. Five-day forecasts were generated by the Holt-Winter method and compared with the observed data. The method was determined to be capable of predicting ionospheric delays and thus enhancing positional accuracy up to 50%. Since the study of the ionospheric delay and forecasting is a new subject in the equatorial region, particularly in Malaysia, a precious model has become essential. Therefore, the study of ionospheric error and forecasting is very important and could help in representing the model error of the data assimilation, which is useful for improving ionospheric models, it is also important to identify and select a suitable model in order to provide GPS data correction to GPS users and to improve the accuracy performance of GPS positioning as well as to assist the future development of satellite-based augmentation system (SBAS) in Malaysia.
In this research, the Holt-Winter method is modified and implemented to forecast ionospheric delay accurately.
The variations in the diurnal, monthly and seasonal ionospheric delay were estimated from the GPS Ionospheric Scintillation and TEC Monitor (GISTM) receiver and compared with the ionospheric delay forecast using the statistical Holt-Winter method over the UKM and Langkawi stations. Besides the differences in ionospheric delay observed during geomagnetic storm disturbances were investigated.

Data and method of analysis
The GPS-TEC data utilized in this work were inferred from the dual-frequency GISTM receivers installed at UKM station (geographic coordinates 2.55°N, 101.46°E and geomagnetic coordinates 2.55°N, 101.46°E), and Langkawi station (geographic coordinates 6.19°N, 99.51°E and geomagnetic coordinates 3.39°S, 172.42°E) in 2011. The GISTM receiver is capable of tracking up to 11 GPS signals at both L1 (1575.42 MHz) and L2 (1227.6 MHz) frequencies. A 50 Hz measurement was used at the amplitude and phase, and a 1 Hz measurement was used for the divergence of satellite's code/carrier. The GPS data was analysed using an elevation angle greater than 20° because the GPS signal suffers from significant fluctuations due to the multipath effect [da Silva and Camargo, 2010], this multipath effect appears with a small elevation angle. The lock time shows the amount of time the receiver was restricted to the carrier phase on L1 and L2 (L1, L2 lock time), the measurement data for L1 and L2 are restricted for over 240 seconds Nouf Abd Elmunim et al.
(4 minutes). This time is required by the GISTM receiver for the detrending High Pass Filter (HPF) to reinitialise when the lock to the carrier phase is lost. The application of a dual-frequency pseudorange and carrier phase allows us the slant TEC (STEC) through the line of sight. The GISTM receiver provided the STEC data every 60 seconds, the STEC was determined from the Receiver Independent Exchange (RINEX) format observation files obtained from the GISTM. The GPS RINEX files were processed using the analysis software for GPS-TEC developed by Gopi Seemala at Boston College's Institute for Scientific Research [Seemala and Valladares, 2011]. The software is available online at http://seemala.blogspot.com. The GPS-TEC analysis software uses the phase and code measurements for both L1 and L2 frequencies in order to calculate the relative STEC values. Then the absolute values of STEC are obtained by including the differential satellite biases published by the University of Bern and the receiver bias that is calculated by minimizing the TEC variability between 02:00 and 06:00 LT [Doherty et al., 2004;Habarulema et al., 2013]. The bias free TEC is given in ASCII files together with the time, elevation angle, azimuth angle, IPP-longitudes and IPPlatitudes parameters. The VTEC is computed by averaging the TEC for individual satellites in view. The converted VTEC in meter unit is termed as vertical ionospheric delays [Abe et al., 2017]. Additionally, the ionospheric delay was approximated using Eq. (1) as shown below (1) where , is the ionospheric delay in metres, with the carrier frequencies f 1 = L 1 and f 2 = L 2 .
The Holt-Winter method is a statistical short-term method that utilises mathematically recursive functions to forecast ionospheric delay trend behaviour [Makridakis et al., 1982]. It uses a time-series with repeated trend and a seasonal pattern to make the forecasting on the assumption that the future will follow a similar pattern. The compared with the Additive model [Elmunim et al., 2015]. Therefore, the Multiplicative model was used in this research. In order to improve the accuracy of the method, the Holt-Winter method equations were modified by identifing some certain smooth coefficients values and added them to the equations. The smooth coefficients founded were 0.9, 0.1, and 0.1 for , and , respectively. The method was modified and tested in UKM and Langkawi stations, resulting in more accurate forecasts than the basic equation results. The Holt-Winter (Multiplicative) method equations are as follows: (2) (4) = ( -1 + -1 ) -1 (5) where is the level; is the trend; is seasonal; is the VTEC, while t is the time period for the , , and components; is the forecasting value of a period ahead; + , is the forecasting period. The smoothing coefficients (0.9), (0.1), and (0.1) are the level, trend and seasonal, respectively; m is the forecast period and s is the seasonal duration. The seasonal initial value and the levels of the seasonal duration are determined using Eqs. (7) and (8) below.
where S 1 is the initial value of the seasonal component; Y 1 and L 1 , are the initial values of the VTEC and the level, respectively; L s is the level of seasonal duration s.
In order to make the forecast, the data of three days were used to forecast the following day. The Mean Absolute Percentage Error (MAPE) was utilised to define the error between the forecast results and the actual data in order to investigate the accuracy and validate the effectiveness of the method using the formula in Eq.
(9) below: where n is the total observed value; PE t the percentage error; F t the forecast value for period t; Y t represents the actual value. Considering that MAPE values less than 10% indicate that the method gives a highly accurate forecast, a range between 10-20% indicates that the method gives a good forecast, a range between 20-50% indicates a reasonable forecast, and lastly a range greater than 50% indicates that the method gives an inaccurate forecast [Lewis, 1982].

Results and discussion
In the equatorial regions, the diurnal variations on quiet days normally vary according to photoionization production and the recombination losses associated with local solar radiation. The diurnal variation of the ionospheric delay at UKM and Langkawi stations during a typical quiet day on 12 May 2011, where the Kp ≤ 1, has been plotted to illustrate the comparison between actual delay (GPS-TEC) and forecast (Holt-Winter) at UKM and Langkawi stations and is shown in Figure 1. The ionospheric delay is displayed on the vertical axis in TECU unit, and local time (LT) is displayed on the horizontal axis in hours. Malaysia's local time is eight hours ahead of universal time (UT). The diurnal pattern of the ionospheric delay observes an obvious single peak. It exhibits a steady increase from morning to a post-noon maximum between 12:00-17:00 LT and then drops to reach a minimum during sunrise at 05:00 LT. Adewale et al. [2011] investigated the diurnal variation of the TEC over Lagos stations in Nigeria. They observed peaks from 15:00 to 17:00 LT, and the minimum at 06:00 LT. The diurnal variation at Bangkok, Thailand, reported by Arunpold et al. [2014], was observed with the trending peak occurring in the post-noon time period between 14:00 and 18:00 LT during 2011. Therefore, this indicates that the peaks of ionospheric delay are associated with the smaller chemical losses at the higher altitudes and the production of solar radiation during the post-noon time [Kelly et al., 1996]. The diurnal behaviour of the actual and forecast ionospheric delays at UKM and Langkawi stations generally demonstrate a good agreement. The Holt-Winter method underestimates the agreement between the actual ionospheric delay during the early morning time. The ionospheric delay at UKM station was slightly higher than Langkawi station. The maximum ionospheric delay at UKM station reached up to 41.9 TECU while at Langkawi station it reached 40.9 TECU. The lower value of the ionospheric delay was 4.7 TECU at both UKM and Langkawi stations.

Monthly variation
The hourly monthly mean of ionospheric delay variations at UKM and Langkawi stations are illustrated in Figure   2 and Figure 3, respectively. In general, it can be observed that the overall actual and forecast ionospheric delays at both stations are comparable in monthly patterns during the year. The ionospheric delays begin to steadily increase from February to reach a maximum during March and April, and then slightly decrease until August before increasing again in September, October, and November. The Holt-Winter method shows a good agreement in the monthly variations with the actual GPS-TEC. The forecast trends were nearly comparable to the GPS-TEC at both the Langkawi and UKM stations, except during the peak hours in the post-noon time, which showed a slight underestimation in March, September, October, and November at UKM station and in March, July, August, and October at Langkawi station. Compared to the Langkawi station, the UKM station recorded higher ionospheric delays during March, August, September, and October. Considering the results obtained, it can be observed that the highest ionospheric delays occurred in the equinoctial months, which are March, April, September, October, and November, whereas the lowest ionospheric delays occurred in the solstice months of January, June, and July. Similar findings over India, Thailand, Kenya, Uganda, and Ethiopia were stated by various researchers [Bhuyan and Borah, 2007;Kenpankho et al., 2011;Olwendo et al., 2012;Arunpold et al., 2014;Tariku, 2015a;Tariku, 2015b;Panda et al., 2015;Kumar et al., 2016].
They reported that the maximum monthly variation was in March and the minimum was during June and July. The highest monthly diurnal peak at both Langkawi and UKM stations was 69.5 TECU, which was likely influenced by geomagnetic disturbances in October, whereas the comparable peak during quiet periods was 32.3 TECU.

5
Characterization of ionospheric delay and forecasting   In the year 2011, a number of moderate geomagnetic storms occurred. August and September were characterised by a number of medium storms. Sub-major geomagnetic storms occurred on 6 August and 26 September, with minimum values of Dst around -113 nT and -103 nT, respectively. These moderate events did not significantly affect the diurnal GPS-TEC and therefore the forecast trends. However, the geomagnetic storm that occurred on 25 to 26 October had a major impact, the details of which will be described in later sections.
In order to validate the forecasting method, the MAPE is used to measure the forecasting errors and determine the forecasting method performance. The diurnal hourly average of MAPE variations, as well as the monthly MAPE averages at UKM and Langkawi stations are shown in Figure 4 and Figure 5, respectively. The maximum MAPE for all the months of 2011 was observed during morning time between 00:00 LT and 07:00 LT, for both stations where the MAPE value was around 8%. In Figure 4 (b) it can be seen that in UKM station the MAPE's maximum monthly average was observed in January, March, and September with a value of about 5.8%, while the minimum was in April and December at around 2.2%. At Langkawi station, the MAPE's maximum monthly average was observed in January, March, and July with a value of about 6%, while the minimum was observed in April and December with a value of about 1.7% as shown in Figure 5

Seasonal variation
To investigate the seasonal variation of ionospheric delay, the diurnal variation of seasonal mean values were calculated for all days of the season, where the seasons are described as summer (May, June, July, August), winter (January, February, November, December), and equinox (March, April, September, October). The seasonal actual and forecast ionospheric delay variations at UKM and Langkawi stations are given in Figure 6. It can be depicted that the maximum value of ionospheric delay was observed during the equinox, where this is due to the sun shines directly over the equatorial region during the equinoctial months causing a strongest ionization over the ionosphere.
The minimum ionospheric delay value was observed during the summer, and that is could be due to unequal heating of the two hemispheres. It is produced by the transport of the neutral constituents from summer to winter hemisphere (hot to cold), which reduces the recombination rate in winter compared to the summer and leads to a relatively increase of electron concentration in the winter than the summer. Moreover, the change of the direction of the neutral wind is another possible cause for this seasonal phenomenon. A meridional component of neutral wind apparently blows from the summer to the winter hemisphere that can reduce the ionization crest value during the summer solstice where it blows in an opposite direction to the plasma diffusion process originating from the Nouf Abd Elmunim et al.
8 Figure 5. The diurnal hourly average of MAPE variation (a) and the monthly MAPE average (b) at Langkawi station. magnetic equator. Consequently, during the equinox meridional wind blows from the equator to polar regions causing a high ionization crest value. Other researchers like Galav et al. [2010], Asmare et al. [2014], and Kumar et al. [2016] studied the seasonal variation in the Indian region, and they found that the maximum seasonal variation occurred during the equinox while the minimum value was in summer. In this research, the maximum seasonal ionospheric delay was observed in equinox at approximately 54.2 TECU at UKM station and 47.6 TECU at the Langkawi station and dropped down to a minimum of 4.7 TECU at both stations. In summer, the maximum ionospheric delay reached up to 35.2 TECU at UKM station and 33.3 TECU at Langkawi station, while it reached a minimum of 1.9 TECU at both stations. Generally, the Holt-Winter method exhibited good agreement with actual GPS-TEC for all seasons and in both stations, except for a slight underestimation around 1.9 TECU which was observed in the peak hours of the summer season at Langkawi station.

9
Characterization of ionospheric delay and forecasting

Disturbed period
The geomagnetic storm occurred between 25 and 26 October 2011 (24 to 25 October in UT) was classified as an intense or major storm, were the Dst value reveals a sharply decrease to -136 nT and the Kp index reachs its maximum to 7 approximately. The diurnal variation of the actual ionospheric delay during the disturbed days from 23 to 28 October 2011 at UKM and Langkawi stations are shown in Figure 8. The geomagnetic storm commenced on 25 October and significantly increased the ionospheric delay from 15:00 LT until the next day at 04:00 LT.
The peak in diurnal variation of the actual GPS-TEC was greatly enhanced on 25 October, when it reached up to 86.6 TECU at UKM station and 80.9 TECU at Langkawi station. Comparing the diurnal ionospheric delay during the disturbed periods, even before and after the geomagnetic storm, with the quiet day of 12 May as discussed earlier in Figure 1, it was found that the ionospheric delay was less during the quiet day than the disturbed periods. During the disturbance of the geomagnetic storm, there is an input of energy into the polar ionosphere, which changes several thermospheric parameters, such as the composition, temperature, and circulation. When the composition changes it directly affects the electron concentration in the F2 region, and then the TEC, while the circulation spreads the heated gas to lower latitudes. Basu et al. [2007] reported that during major geomagnetic storms the 11 Characterization of ionospheric delay and forecasting currents associated with inner magnetospheric electric field directed in the dusk-to-dawn direction can no longer shield the mid-latitude and equatorial latitudes from the high-latitude electric fields, and this leads to instantaneous penetration of electric fields from high latitude to the middle and the equatorial ionosphere. Consequently, the particle transport and the prompt penetration of the high-latitude electric field to lower latitude which travel equatorward with high velocities during the storm can cause an increase in the TEC value and the ionospheric delay.
The diurnal variation of the actual ionospheric delay (GPS-TEC) during the geomagnetic storm disturbance is shown in Figure 9. continued to increase even when the IMF Bz attained its maximum positive value [Davis et al., 1997]. The Dst started to decrease from 05:00 LT on 25 October (21:00 UT, 24 October) and dropped to a minimum value of -137 nT at 09:00 LT (01:00 UT) on 25 October. After the commencement of the geomagnetic storm, the Dst trend did not recover to its initial value until three days later. At high latitudes, the electric field penetration leads to the formation of ionospheric electron density irregularities in the equatorial region. Therefore, the greatest decrease in the Dst value was observed when the sudden intensification of the ring current occurred [Basu et al., 2007]. The positive value of the Dst mainly occurs due to the compression of the dayside magnetosphere during the initial phase of the geomagnetic storm, and the negative value is due to magnetic reconnection and the formation of the ring current during the main phase of the geomagnetic storm. Several studies have discussed the perturbations during a geomagnetic storm in diurnal variations [Rama Rao et al., 2009;Kumar et al., 2014]. Zou et al. [2013] used a GPS-TEC to observe the same storm on 25 October 2011, and they reported that the ionospheric delay increased at a higher rate compared to that seen during the quiet time. This effect was more clearly observed at lower latitudes than at higher latitudes. Other researchers studied the effects of the same geomagnetic storm on 24-25 October 2011 on the ionospheric TEC for two East African stations; they found that for both stations, the diurnal variations showed a decrease during the storm period [D'ujanga et al., 2013]. October 2011 for both the UKM and Langkawi stations. The vertical axis shows the ionospheric delay in TECU while the horizontal axis shows the hourly mean of the disturbed period days. The forecasting variations in disturbed period ranged from 3.8 -80.9 TECU at the UKM station and from 2.8 -80 TECU at the Langkawi station. In both stations, the maximum value for ionospheric delay forecasting was observed during the geomagnetic storm day, 25 October. The percentage of the MAPE in Figure 10 (b) shows that the UKM station had a comparatively lower forecasting error than the Langkawi station. The maximum forecasting error for the UKM station was 7.3%, while it was 8% at the Langkawi station. The MAPE value was slightly higher during the disturbed period than the quiet period due to the increase of ionospheric delay during the intense geomagnetic storm that occurred on 25-26 October. This increase affected the forecasting and led to a slightly higher forecasting error.
However, based on the result of the ionospheric delay variability during the quiet and disturbed periods, it can be observed that there were no significant changes in the trends during these two periods. The same finding was presented by Habarulema et al. [2013], for their study of the African equatorial region; they observed two severe geomagnetic storm events which occurred in November 2004. Nonetheless, even a slight increase must be taken into consideration, as this has the effect of comparatively increasing the forecasting error compared with the results obtained during the quiet periods. Thus, this will lead to a slight decrease in the accuracy of the forecasting method during the disturbed periods.

Conclusion
This research presents the variability of ionospheric delay and forecast using the statistical Holt-Winter method in Malaysia over UKM and Langkawi stations for the year 2011. The results show that the diurnal variability of ionospheric delay has minimum value during sunrise at 05:00 LT and maximum during the post-noon time mostly between 12:00-17:00 LT. The Holt-Winter method forecast displayed a good agreement, with a slight underestimation compared to the actual GPS-TEC observed in both UKM and Langkawi stations for the entire year.
The maximum and minimum monthly ionospheric delay were observed on March and July, respectively. However, the ionospheric delay had a significant increase in October due to the geomagnetic storm effect. UKM station showed a higher ionospheric delay than Langkawi station. The ionospheric delay was lower during the summer season and higher during the equinox season. The maximum ionospheric delay forecast was observed during the geomagnetic storm day on 25 October in both stations. Comparing the MAPE values during quiet and disturbed periods, it can be observed that the MAPE in disturbed periods is slightly higher than the MAPE during quiet periods. However, since it is still less than 10% it is considered to be an accurate forecasting result. Langkawi station had a relatively higher MAPE than UKM station. The obtained results showed the effectiveness of the Holt-Winter method in forecasting ionospheric delay, providing good results in the diurnal, monthly and seasonal variations, and also during the disturbed period. These results could be useful in representing the model error using data assimilation.
Consequently, this could be also useful for developing the regional ionospheric delay model for better understanding and forecasting over Malaysia. Besides it might help to improve the accuracy and the performance of GPS positioning in equatorial regions, and particularly in Malaysia.