Single-Frequency Time-Step Ionospheric Delay Gradient Estimation at Low-Latitude Stations

The irregularity of the local-area ionospheric delay is a primary impediment for Ground-Based Augmentation System (GBAS) services. Excessive ionospheric delay gradients may degrade aircraft positioning for high precision landing systems. Therefore, the spatial gradients of the nominal background ionosphere must be studied as their statistics will be sent to the approaching aircraft. For the well-known station-pair method, ionospheric delay gradient estimation requires at least 2 Global Navigation Satellite System (GNSS) reference stations. This method can be applied to both single or dual-frequency GNSS receivers. However, when the GNSS stations are far apart, it is not suitable for estimating the ionospheric delay gradients at short baselines, and the time-step method is an attractive alternative. In this work, we propose a single-frequency time-step method for ionospheric delay gradient estimation. Careful baseline length selection is needed, due to ionospheric piercing point movements. We applied our method to GNSS data in 2014, at the peak of the 24th solar cycle, and showed that the standard deviations of the vertical ionospheric delay gradients were comparable to those derived from the dual-frequency time-step method. The standard deviations of vertical ionospheric gradients, $\sigma _{\mathrm {VIG}}$ , ranged between 4 and 6 mm/km. The $\sigma _{\mathrm {VIG}}$ values around the equinoxes were ~1.5 mm/km greater than at other times.


I. INTRODUCTION
The Global Navigation Satellite System (GNSS) is commonly used in the aircraft navigation system, during many phases of a flight. The International Civil Aviation Organization (ICAO) has proposed a Ground-Based Augmentation System (GBAS) [1] as an ICAO standard for the auto-landing system in an aircraft, on or after the approaching phase, within 42.6 km (23 nautical miles) from an airport [2]. In this system, the aircraft position accuracy can be improved by using the differential correction from ground reference stations at the airport. Three or four multi-frequency GNSS receivers must be installed around the area for the GBAS service. However, variation in the local area ionospheric delay is a primary impediment for the GNSS signal quality [3]. The aircraft positioning errors, both vertical and horizontal, The associate editor coordinating the review of this manuscript and approving it for publication was Poki Chen . estimated with the ground correction service, have to be maintained in an acceptable range [4], so performance monitoring of the positioning systems is needed.
In GBAS, statistics of the ionospheric variations [3] is one of the parameters in the performance monitoring systems. The understanding of normal vertical ionospheric delay gradients (σ VIG ) is important for GBAS operation. The parameter σ VIG is used to calculate the protection level, which bounds the error in the usual fault-free conditions. Therefore, it ensures safety in common operating condition. Moreover, σ VIG has sometimes been used to prevent unsafe operation in disturbed conditions, by inflating it in a GBAS system [5]. Previously, ionospheric parameters were studied in both quiet and disturbed events for GBAS operation [6]- [9]. During the quiet condition, in the low-latitude region (3 countries), the ionospheric delay gradients in the North-South direction were higher than the East-West direction, due to the Equatorial Ionospheric Anomaly (EIA) [10]. From the disturbed ionospheric studies in the Brazilian region [9], high ionospheric delay gradients were observed within ±35 degrees of the magnetic latitude, mainly caused by the nocturnal Equatorial Plasma Bubbles (EPBs). In March 2012, 2014, and 2015, the ionospheric delay gradients in reached 505.2 mm/km in Singapore [11], 850.7 mm/km in Brazil [12], and 308 mm/km in India [11]. These gradients were significantly larger than the normal condition and showed how vital the ionospheric monitoring for GBAS operation in low-latitude regions is. Even in quiet times, statistics of the local delay gradients must be transmitted from runways to the approaching aircraft. Local monitoring with accurate statistics helps increase the availability of the GBAS system. Current ionospheric delay gradient estimation methods, for example, station-pair, time-step, and the mixed-pair methods, are shown in Table 1. In the station-pair method, the delay gradients are estimated from the total electron contents, TECs, obtained from the two frequencies, used in GPS receivers, which always use an L1 (1575.42 MHz) and, sometimes, an L2 (1227.60 MHz) signal, between two stations, but, in reality, we may not have the L2 signal from the Global positioning system (GPS) receiver, due to ionospheric irregularity constraint of L1 only GPS receivers on aircraft. In [13], [14], the single-frequency station-pair delay gradient estimation was proposed based on the L1 frequency with the Kalman filter and the LAMBDA method. The satellite elimination [15] was also applied to this method to increase the success of the estimation. However, for most areas in Thailand, the GNSS receivers are widely separated (more than 20 km baselines) since they are not originally intended for delay gradient monitoring. The miss detection of the high delay gradients over short baselines remains. However, the time-step method [16], [23], in which the delay gradient estimated from a single station at different time instants, requires dual-frequency observations to estimate the delay gradients. However, the time-step method cannot be directly applied to single-frequency GNSS receivers due to the baseline variation issue. Hence, in this article, a single-frequency time-step method is described. The time index management for IPP-pair selection was applied. The baselines at the ionospheric heights were appropriately adjusted for delay gradient estimation at each epoch.
Our method was used at two stations near Suvarnabhumi International Airport, Thailand and validated by comparing with single-frequency and dual-frequency station-pair methods. This article is organized as follows. The existing methods for delay gradient estimation and our new method are described in Section 2. The experimental setup and the parameters of delay gradient comparison are described in Section 3. The compared results are shown in Section 4. Finally, we conclude in Section 5.

II. IONOSPHERIC DELAY GRADIENT ESTIMATION A. EXISTING METHODS
A point, the ionospheric pierced point, IPP, is defined as the point at which the signal passes into the ionosphere. The delay gradients refer to the difference between two slant delays, divided by the separations of the IPPs. The key to achieve this is receiving the satellite signals from the two different locations or at different time instants. Three methods are used, i.e., the station-pair [24], the time-step [19] and the mixed-pair methods [20], shown in Fig. 1.  [24], (b) time-step [19] and (c) mixed-pair [20].
In Fig. 1, for the station-pair method [24], two GNSS stations (R1 and R2) are installed at two permanent locations. The GNSS signals from both receivers are used to obtain the delay gradients. Since the GNSS receivers are permanent, the baseline length, D ground , at ground level is constant. Next, for the time-step method [19], only one GNSS receiver is used. The GNSS signals from each satellite are received at two-time instants, T 1 and T 2 , for any separation of interest. If the time difference, T = T 2 -T 1 , corresponds to a short baseline D, between 1 and 40 km, the estimated delay gradients can be used for the GBAS applications [16]. Lastly, for the mixed-pair method [20], IPP pairs selected from any satellite and receiver with various baseline lengths, are used to obtain the delay gradients. However, this method requires two distinct frequencies, and the varying IPP pairs prevent the VOLUME 8, 2020 single-frequency approach to fix the integer ambiguities [25]. Although the IPP pairs, from the mixed-pair method, include more combinations than the station-pair method, in local area ionospheric monitoring, the station-pair method is sufficient.
Both station-pair and time-step methods estimate the delay gradients from the signals received from one satellite at a time. The GNSS satellite signals to the receivers intersect at different IPP locations. For example, in the station-pair method, two IPPs are at points 1 and 2. However, the main difference of delay gradient estimation is the number of the GNSS receivers used in each method. If there is only one GNSS receiver in the area, the time-step method can be used to estimate the delay gradients.

1) DUAL-FREQUENCY IONOSPHERIC DELAY GRADIENT ESTIMATION USING THE STATION-PAIR METHOD
The ionospheric delay gradients, ∇I g (t), corresponding to a signal path from a satellite, g, can be estimated from the difference between the total electron content from two stations divided by the baseline length, D ground (km) [24], i.e., where f is the frequency of the GNSS signal (i.e., 1575.42 MHz for the L1 GPS signal), STEC g R1 (t) and STEC g R2 (t) are the total electron content along the slant paths to stations, R1 and R2.
Typically, the delay gradients based on dual-frequency GNSS data, in (1), can be estimated using either the stationpair or time-step methods. However, in the receivers at the ground station, occasionally, the L2 signal has relatively lower signal-to-noise ratios (SNR), than the L1 counterpart. Additionally, aircraft mostly use an L1-only receiver, so the second frequency is unavailable. Therefore, the single-frequency gradient estimation is an alternative approach under low SNR situations.

2) SINGLE-FREQUENCY GRADIENT ESTIMATION USING THE STATION-PAIR METHOD
For the single-frequency gradient estimation [13], two types of pseudo-range measurements are used to measure the distance from a receiver to a satellite. The code measurement, ρ, and the carrier phase measurement, , of a single station can be generally expressed as and where r is the true distance between the receiver and the satellite in meters (m), the superscript g represents the pseudorandom number of the satellite, b is the receiver clock bias (m), B is the satellite clock bias (m), δI is the ionospheric delay (m), T is the tropospheric delay (m), λ is the GNSS signal wavelength, N is the integer ambiguity (cycles), and ε and ε ρ are the measurement noises, relative to the multipath effects, from the code and the carrier phase measurements. Firstly, to remove the ionosphere term included in the pseudo-range and the carrier phase, the ionosphere-free combination,L, is computed from The single difference (SD) is defined as the difference of code or carrier phases between stations R1 and R2,L g SD , which helps remove the satellite clock bias terms, can be computed fromL Therefore, from (2) to (5), we havẽ Note that b SD and N g SD account for the differences, in receiver clock bias and integer ambiguity of the two paths. Similarly, by using the SD combination described above, the single difference of geometry-free carrier-phase (˜ g SD ) can be expressed as In (5) to (7), for short baselines between receivers [26], the tropospheric delay is assumed to be constant, within a large area. Assuming that the single difference of tropospheric delay, T SD = 0, (5) and (7) together can be expressed in a matrix, andḃ SD is receiver clock bias drift, : n × 1, n is the number of the observed satellites, I is an n × n identity matrix and O is n × n zero matrix. From (8), the vector, δI SD , the desired delay gradients from two GPS receivers and other unknowns, can be estimated by Kalman filtering [13] and fixed the ambiguity by LAMBDA method [25]. The baseline computed between two receivers is used to obtain the vertical delay gradients. Then, the standard deviation is computed from one-day delay gradient data. However, an actual distribution, estimated from the observation data, has non-Gaussian tails. In the GBAS standard [1], the zero-mean Gaussian distribution is used as an error model in the protection level computation. A nominal sigma, one σ , of a zero-mean Gaussian distribution, must be inflated to make a system safe, i.e., to reduce the probability of deviations to < 10 −7 . This inflated sigma of vertical delay gradient is called σ VIG [8]. From the reference method [13], the baseline was computed at the ground level (Earth's surface), as shown in Fig. 2. This baseline at ground level does not directly represent the ionospheric condition. To improve the gradient estimation, this baseline must be computed from two IPP locations. Hence, here, the baseline computed at both heights was used to analyze its effect on the gradient estimation.
The single-frequency station-pair method requires at least two stations. In this article, we describe a single-frequency time-step method to estimate the delay gradients based on a single GNSS receiver.

B. OUR SINGLE-FREQUENCY TIME-STEP METHOD GRADIENT ESTIMATION
The previous section discussed existing gradient estimation techniques, including the station-pair and the time-step methods. However, the single-frequency time-step method ionospheric delay gradient has not been implemented yet. The baseline variation was reduced before the Kalman filter was applied to the estimation. These baselines in the time-step method are described in Fig. 3.
In the time-step method, the baselines at times, T 1 and T 2 , corresponding to the two IPPs must be computed. Normally, the IPP coordinates can be obtained from the satellite position and the receiver position in both the cartesian (s ECEF , r ECEF ) and the geographic coordinate (r lat , r lon ) system. The varying baselines were used to compute the delay gradients.
The time-step method assumes that the ionosphere is stationary within the time interval. If the ionosphere changes with a speed similar to the IPP speed, the observations needed to calculate the high delay gradient may be missed. Therefore, in this work, the single-frequency time-step method was used to compute the background ionospheric gradient statistics. These background values do not change rapidly, compared with the speed of IPP movement.

1) IONOSPHERIC PIERCE POINT (IPP) COORDINATES
At any given satellite and receiver position, the IPP coordinate in a geographic coordinate system can be obtained from [27], i.e., and where r lat is the latitude of the receiver, R e is the earth's radius, h is the ionospheric height, Az and El are the azimuth and the elevation angles from a receiver to a satellite. The latter two can be obtained from and The transformation matrix RL = RL x RL y RL z can be obtained from where and, (16) as shown at the bottom of the next page, where s ECEF and r ECEF are the vectors of satellite and receiver positions in the cartesian coordinates. VOLUME 8, 2020 The IPP positions from a given satellite and receiver position can be used for the baseline estimation. However, the original haversine formula computes the ground baseline. When used in the time-step method, the baseline from the ground and the IPP height should be considered. The next section describes the obtained baseline at both levels.

2) ADJUSTED IPP DISTANCE AT THE IONOSPHERIC HEIGHT
Normally, at each instant, the baseline from two different IPPs is shown in Fig.3. From the haversine formula [28], the baseline at the ground level from these IPP locations can be calculated from where R e is the earth's radius (6721 km) and The parameter a can be computed from when ϕ is the difference of latitudes at two locations, λ is the difference of longitudes at two locations, ϕ 1 and ϕ 2 are the longitudes at both IPP locations. However, the baseline must be computed at the IPP level. Therefore, by using the haversine formula at the height of 350km (lower level of the ionosphere), the adjusted IPP distance (d u ) can be computed

3) SELECTION OF TWO IPPS CLOSE TO THE DESIRED BASELINE
From Fig. 4, in the time-step method, if a GNSS station received the signal at time, T 1 , that passed through the IPP coordinate, IPP lat,lon,1 , with range measurements of the code, ρ T , and carrier-phase, φ T 1 , and, at time T 2 , the IPP moved to IPP lat,lon,2 , with corresponding range ρ T 2 and code φ T 2 ,, then the time difference, T = T 2 − T 1 , must be smaller than the temporal change of the ionosphere (about 10 min.), the delay gradients can be obtained from [16]. The benefit of the time-step method is that the delay gradients can be obtained from only one dual-frequency GNSS station, hence, it is suitable for a baseline experiment with a single receiver in an area of interest. Moreover, the baseline variation, due to the satellite speed, does not affect the dualfrequency estimation, due to the time-independent estimation. If a specific baseline is required, baseline selection can be applied after using dual-frequency estimation. However, to apply the time-step method to single-frequency estimation, the SD of ionospheric delay in the previous and the predicted states in the Kalman filter should be selected from the same baseline length to reduce estimation error. Hence, a method to adjust T for obtaining the desired baseline is described as shown in Fig. 5.
The flowchart in Fig. 5 can be used to determine the time difference that provides a baseline close to the threshold. Firstly, if the signal from the satellite is observed at time, T k , the preliminary baseline is computed from two IPPs at times, T k , and at the previous epoch, T k−1 . Additionally, by using the adjusted Haversine formula, (20), the baseline is computed at the IPP height. Next, if the initial baseline is less than the specified threshold, the time difference will be extended until it can find suitable epoch, T k−s , that leads to a baseline greater than the threshold. Then, the pseudo-range at T k and T k−s is applied for the single-frequency delay gradient estimation.

4) MODIFIED SINGLE-FREQUENCY DELAY GRADIENTS ESTIMATED WITH THE TIME-STEP METHOD
From the new time selection technique, a suitable baseline was found by using pseudo-ranges from two proper time indices. In the time-step method, since we have a single GNSS station, the SD combination must be modified to compute the difference in delays at times, t − t and t instead. For the carrier-phase, the SD in (4) and (6) becomes 201520 VOLUME 8, 2020 and, similarly for the code delays, where t is the time difference that gives the baseline, d, between two IPPs, close to the threshold. The modified SD combinations from (21) to (23) can be used to compute the delay gradients. based on the time-step method (single station). Note that at each epoch, due to satellite movement, the true range, r, must be re-estimated.

III. EXPERIMENTAL SETUP
To compare the performances of the ionospheric delay gradients from both methods, we used two GNSS stations (KMITL at N 13.7276, E 100.7724 and STFD at N 13.7356, E 100.6611, with KMITL bearing 105.08 • from STFD) in the testing area. Both GNSS antennas were located on building rooves, as shown in Figs. 6 and 7. These buildings are the tallest buildings within a radius of 10 km. The multipath is minimized, due to the antenna heights, no surrounding buildings can block signals. The calibrated positions of both antennas were obtained from the precise point positioning  (PPP) technique [29], [30]. The KMITL location was used for both the dual-frequency time-step method and our singlefrequency time-step method.
For the baseline estimation in the time-step method, the ground baseline from the STFD to KMIT stations was used as the reference. To select the baselines, at the IPP height of about 12 km, GPS satellite data on day 103 in 2014 were used, due to the quiet ionospheric conditions. The GPS data of 32 satellite indices or PRNs from 00:00:00 UTC to 23:59:59, 86400 epochs, was used in the baseline estimation.
To provide the baselines between two IPPs from the timestep method, the baseline from the station-pair method was selected as a threshold. The KMITL and the STFD stations were the reference baseline -about 12 km. As a final point, this distance was used to decide the time differences applied to each satellite in the time-step method for a similar baseline.
An example of GPS data on DOY 103 in 2014 was selected for the gradient estimation and the baseline estimation by the station-pair and our new single-frequency time-step methods. The ionospheric condition on this day was verified as quiet based on the Rate of TEC change Index (ROTI) [31]. Finally, all days of 2014, at the peak of the 24th solar cycle [32], are used to compare the σ VIG from various methods.

IV. RESULTS AND DISCUSSIONS A. ADJUSTED IPP DISTANCES IN THE TIME-STEP METHOD
Baseline lengths computed from the time-step method shows how the baseline length can vary in one day. Normally, most VOLUME 8, 2020 of the baseline lengths for each ionospheric analysis can be made by changing to sampling over a full day. Moreover, at each sampling time, the performance of the baseline length selection can be measured from the baseline length variation. This value shows how much data must be removed before using the ionospheric delay gradient estimation. To see the performances, when the baseline length selection was not applied, the results from each sampling time are shown in Fig. 8. Fig. 8 shows the baseline results from 32 GPS satellites plotted in different colors. The x-axis is the UTC (hour) and the y-axis is the baseline length in the time-step method (km.). The adjustment of the sampling times from the reference method can vary the lowest baseline in the calculation of delay gradients. For example, the 1-min sampling time leads to the lowest baseline, ∼3 km, as shown in Fig. 8-a, whereas at 5-min sampling time, the lowest baseline is ∼15 km, as shown in Fig. 8-b. However, both sampling times contain unwanted baselines, due to relative motion between the satellites and the receivers. The 1-min sampling leads to an undesirable baseline of 18 km, which increases to 80 km with 5-min sampling. Therefore, to show the improvement in the baseline length variation, the same dataset for the time-step method was used for the IPP-pair selection, as shown in Fig. 9.
From the selection of the baseline length result, the baseline lengths were adjusted from T = T k − T k−s in the time-step method. The threshold of the distance was set to 12 km. Then, the T was adjusted until the baseline was as close to 12 km as possible. Then T to each satellite was used to obtain the delay gradients. Most of the baselines were ∼12 km, with less than a 600 m error. No significant-high baseline length must be removed before using the delay gradient estimation. All of these baselines were optimal for the delay gradient estimation and compared with the reference single difference station-pair method.  Fig. 9, different sampling times in the time-step method caused the baseline to vary. According to the adjusted IPP results, these baselines were computed in the range from 12 to 12.6 km. Therefore, to show the relationship between the ground and the IPP baseline, the one-day satellite positions were used in the analysis, shown in Fig. 10. Fig. 10 shows that the baselines, computed from both ground and IPP height, in the time-step method are linearly related. The azimuth and the elevation angle of the satellite did not affect the computed baselines at both heights. The baselines of 12 and 12.7 km at the ground were adjusted to about 12.66 and 13.27 km at the IPP height, i.e., at the IPP height, they were 5% higher than on the ground. We concluded that the previous delay gradients, estimated from the time-step method, with the ground baseline, should also be adjusted by reducing them by ∼5%, as we change to gradients estimated with the IPP baseline.

C. GROUND BASELINE LENGTH AND THE IONOSPHERIC BASELINE LENGTH IN THE STATION-PAIR METHOD
From Fig. 2, baselines computed from the haversine formula (17), at the ground level formula, D ground = 12.0777 km. However, when the haversine formula was applied at the IPP height, a realistic satellite position is required. The 71.84 • elevation and 350 • azimuth, from the satellite navigation file, was used. The computed baseline from this setup is ∼11.8223 km. When compared to the haversine formula at ground level, the new baseline was 2.11% lower than the reference one. This value was used to increase the delay gradient estimate. Hence, more data should be analyzed. Since the baseline output can vary with the satellite position, and in any direction, the analysis from the station-pair method, with different azimuth and elevation angles, over a day of observation, is shown in Fig. 11.   FIGURE 11. Baseline estimation at the IPP height from the station-pair method from satellite positions over a day. Fig. 11 shows that baselines computed from the West and East direction, below 30 • elevation angles, differ significantly between the ground and IPP levels. The baselines computed at 30 • is ∼15% different. They increased to ∼50% difference, for elevation angles < 15 • . However, the baselines computed from the North and South directions were more consistent; here baseline differences, for elevations from 0 to 30 • , were lower than 8%. Therefore, previous delay gradients computed from the station-pair at the ground location, from satellite positions above 30 • , do not need to be adjusted, for the estimation at the IPP level. However, if delay gradients, estimated from below 30 • satellite elevations, are to be obtained, results should be adjusted by the difference of baseline estimation from the ground and the IPP levels indices.

D. IONOSPHERIC DELAY GRADIENTS FROM THE REFERENCE SINGLE-FREQUENCY STATION-PAIR METHOD
This part describes the delay gradients estimation using the reference single-frequency station-pair method, shown in Fig. 12. The y-axis shows the ionospheric delay gradients (mm/km). Each color represents the delay gradient estimated from a different satellite. Positive values show that ionospheric delays from the East are greater than from the West. The x-axis shows the UTC time. Unsuccessfully resolved integer ambiguities in the LAMBDA method were removed from the figure. Finally, the delay gradients in this plot are vertical values, converted from the slant gradient.
In Fig. 12., delay gradients lie in the range between −20 to 20 mm/km and represent a quiet time. To compute σ VIG , the actual standard deviation was inflated by a factor of 1.05. The σ VIG from this method was about 5.42 mm/km. The σ VIG estimated from the single-frequency method showed that these delay gradients represent quiet ionospheric conditions. The delay gradients from the station-pair method showed East-West symmetry all day, shown by an average close to 0 mm/km. These parameters were used to compare with our new time-step method estimation.

E. DELAY GRADIENTS ESTIMATED FROM OUR NEW SINGLE-FREQUENCY TIME-STEP METHOD
For the time-step method, if the delay gradient was estimated at time T k , the IPP pair was selected from times T k and T k−s . VOLUME 8, 2020 The IPP direction can be measured from IPP location at time T k−s to T k . However, the directions, computed from each IPP pair, cannot be assigned initially, they are computed only after selecting the IPP pairs. If analysis of the delay gradient in a specific direction is desired, it can be chosen only in the last procedure. To show the direction of each pair, in the Thai region, a histogram of one-day IPP directions is shown in Fig. 13. The x-axis shows the direction of the IPP pair from 0 to 359 • . The y-axis shows the frequency in each direction.  The IPP pairs computed from a single day observation were not in uniformly distributed directions, but mostly lie in the North and South directions, i.e., frequencies within 10 • of the North-South axis, in IPP direction, constituted about 10 5 epochs, but less than 1% of all data was in the East-West directions. There was clearly a disproportionate number in the North-South direction. Moreover, if a raw delay gradient estimate were used (unassigned direction), a single direction of the delay gradients would not be observed. In Fig 14, South-to-North directions are indicated with -ve signs. Using our new time-step method, delay gradients range from −10 to 30 mm/km. Delay from the South (low latitudes) were higher than from the North (high latitudes). This statistic indicated that the delay gradients from high latitudes are mostly lower than in equatorial latitudes. The standard deviation of delay gradients from the single-frequency time-step method was ∼4.89 mm/km, comparable that derived by the previous single-frequency station-pair method. Moreover, to show the correlation of the delay gradients estimated by the new and previous methods, the East-West direction (with outliers eliminated) of the delay gradients from the two methods are shown in Fig. 15. The East-West direction from our method was considered, since IPP pairs from the previous station-pair method were obtained from this direction only.
Since the IPP pair directions from the time-step method mostly lie in the North-South direction, fewer than 100 samples of the preferred East-West direction were left for comparison. From this East-West data, the correlation coefficient, r, between the single-frequency station-pair method and our new single-frequency time-step method is ∼0. 25. This shows that the delay gradients, estimated from the single-frequency, in both station-pair and the time-step methods, exhibited the same trend for the same direction of delay gradients. Finally, σ VIG estimated from these methods in 2014 are shown in Fig. 16. Correlations between the new method and the reference methods are shown in Table 2. To fairly compare the time-step and other methods, since the station-pair method provides only East-West delay gradients, the angles of IPP between 45 • and 135 • degrees were compared. The expected σ VIG is about 5 mm/km. From Fig. 16, in 2014, the new single-frequency time-step method yielded gradients that were 30% lower than the reference single-frequency station-pair method. Since the data at each time was from the same GPS receiver, the effects of the receiver biases were almost completely canceled. Therefore, the new single-frequency time-step method provided σ VIG similar to the dual-frequency time-step method, especially during the March and September equinoxes. The estimated σ VIG during both equinoxes was about 1.5 mm/km greater than other seasons. During the June and December solstices, both dual-frequency and our single-frequency time-step method yielded the lowest σ VIG . This was consistent with our previous measurements in the seasonal analysis [21]. Table 1 shows the correlation coefficients between the new method and the dual-frequency time-step method (r = 0.68) and the single-frequency station-pair method (r = 0.32). Thus, our method provided σ VIG similar to the dual-frequency time-step method. However, the correlation coefficient between the new method and the single-frequency station-pair method was lower. We attribute the choice of only East-West IPP pairs in the new time-step method and their coordinates were mostly different from the IPP pairs in the single-frequency station-pair method and the ionospheric conditions may have been slightly different.
Finally, since the time-step method required only one GNSS station to produce the delay gradients, ionospheric monitoring from the time-step method is more available than with the reference station-pair method. The existence of the other nearby stations can provide redundancy to ensure that the delay gradients are accurate. Our method is also useful on single-frequency mobile phone, Ublox GNSS receiver [33].

V. CONCLUSION
We described a single-frequency time-step method to estimate the ionospheric delay gradients which needed only a single GNSS station. We tested it around Suvarnabhumi airport, Thailand. With the new single-frequency time-step method, most baselines were not over 12.6 km. The baseline variation improved to below 600 m within the 1-min data sampling time. These baselines were optimal for comparison with the delay gradients estimated from the single difference station-pair method. Next, when the stations from the station-pair method were oriented East-West, the baseline estimated at the ground and IPP levels differed for elevation angles < 30 • : baseline differences were about 15% different at that angle. and they increased to about 50% different for elevation angles < 15 • . However, in the time-step method, the azimuth and the elevation angle of the satellite did not affect the estimated baselines at both heights. The baselines from both ground and IPP heights were linearly related for any satellite positions. Delay gradients from higher latitudes were mostly less than at equatorial latitudes. The single-frequency time-step method also yielded lower delay gradients ∼30% than the reference single-frequency station-pair method. The dual-frequency and our new single-frequency time-step method estimated σ VIG around the equinoxes to be 1.5 mm/km greater than at other times. However, IPP pairs computed from one-day observations were not equally distributed in all directions. The relative motion between the satellites and the receivers limits the time-step method to be only suitable for the delay gradient monitoring in the North-South direction. Finally, the preliminary results of our new single-frequency time-step method would be useful for a standalone GNSS station to detect disturbed ionospheric conditions. However, to identify these conditions, the single-frequency cycle slip detection must be implemented along with time-step approach. It is an important future work in ionospheric monitoring using single-frequency GNSS receivers.