Generation of Ionospheric Delay in Time Comparison for a Specific GEO Satellite by Using Bernese Software

Abstract

Time comparison is necessary for the verification and synchronization of the clock. Two-way satellite time and frequency (TWSTFT) is a method for time comparison over long distances. This method includes errors such as atmospheric effects, satellite motion, and environmental conditions. Ionospheric delay is one of the significant time comparison error in case of the carrier-phase TWSTFT (TWCP). Global Ionosphere Map (GIM) from Center for Orbit Determination in Europe (CODE) is used to compare with Bernese. Thin shell model of the ionosphere is used for the calculation of the Ionosphere Pierce Point (IPP) between stations and a GEO satellite. Korea Research Institute of Standards and Science (KRISS) and Koganei (KGNI) stations are used, and the analysis is conducted at 29 January 2017. Vertical Total Electron Content (VTEC) which is generated by Bernese at the latitude and longitude of the receiver by processing a Receiver Independent Exchange (RINEX) observation file that is generated from the receiver has demonstrated adequacy by showing similar variation trends with the CODE GIM. Bernese also has showed the capability to produce high resolution IONosphere map EXchange (IONEX) data compared to the CODE GIM. At each station IPP, VTEC difference in two stations showed absolute maximum 3.3 and 2.3 Total Electron Content Unit (TECU) in Bernese and GIM, respectively. The ionospheric delay of the TWCP has showed maximum 5.69 and 2.54 ps from Bernese and CODE GIM, respectively. Bernese could correct up to 6.29 ps in ionospheric delay rather than using CODE GIM. The peak-to-peak value of the ionospheric delay for TWCP in Bernese is about 10 ps, and this has to be eliminated to get high precision TWCP results. The $10^{-16}$ level uncertainty of atomic clock corresponds to 10 ps for 1 day averaging time, so time synchronization performance needs less than 10 ps. Current time synchronization of a satellite and ground station is about 2 ns level, but the smaller required performance, like less than 1 ns, the better. In this perspective, since the ionospheric delay could exceed over 100 ps in a long baseline different from this short baseline case, the elimination of the ionospheric delay is thought to be important for more high precision time synchronization of a satellite and ground station. This paper showed detailed method how to eliminate ionospheric delay for TWCP, and a specific case is applied by using this technique. Anyone could apply this method to establish high precision TWCP capability, and it is possible to use other software such as GIPSYOASIS and GPSTk. This TWCP could be applied in the high precision atomic clocks and used in the ground stations of the future domestic satellite navigation system.

1. INTRODUCTION

Time comparison is necessary for the clock verification and synchronization (Hwang et al. 2016). Since the optical fiber link has restriction in intercontinental case, satellite is used for overseas connections (Fujieda et al. 2014). Two-way Satellite Time and Frequency Transfer (TWSTFT) is a method for time comparison over long baselines (Parker & Zhang 2005). Piester et al. (2007) applied this method for the cesium fountain clock comparison between Asia and Europe. Hwang et al. (2016) analyzed Asia link between Korea Research Institute of Standards and Science (KRISS) and National Institute of Information and Communications Technology (NICT) by using this technique. However, TWSTFT shows instabilities due to error sources such as atmospheric effects, satellite motion, and environmental conditions (Piester et al. 2007). Therefore, these errors have to be eliminated for the high precision time comparison. Among error sources, ionospheric delay is considered significant in the carrier-phase TWSTFT (TWCP) compared to code due to high precision (Fujieda et al. 2014).

It is important for the TWCP to eliminate ionospheric delay in order to get better results. This is because the clock comparison requires high precision. In clock development, the uncertainty of the optical clock is much better than the low level of the cesium fountain clock (Yu et al. 2017). In addition to this, since the BeiDou system has used TWSTFT as the time synchronization method for ground stations (Han et al. 2011), the TWCP could be used in the development of the future domestic satellite navigation system. Moreover, organizations around the world using TWCP could establish more precise time comparison by applying ionospheric delay elimination technique presented in this paper.

Ionospheric delay does not cancelled out due to different frequency between uplink and downlink in TWSTFT (Fujieda et al. 2014). Total Electron Content (TEC) could be estimated on the condition that two downlink frequencies available like Global Positioning System (GPS). However, it is not possible for TWCP that other sources like Global Ionosphere Map (GIM) or regional Vertical Total Electron Content (VTEC) map have to be used. Ionospheric delay is calculated by using global TEC map of the International GNSS Service (IGS) (Parker & Zhang 2005). Piester et al. (2007) used GIM from Center for Orbit Determination in Europe (CODE). Fujieda et al. (2014) have taken into account regional VTEC maps from both the Royal Observatory of Belgium and NICT.

It is necessary to calculate slant total electron content (STEC) in order to get ionospheric delay for uplink and downlink in time domain. First, ionosphere pierce point (IPP) is calculated between stations and Eutelsat 172A, and then VTEC at that point is calculated by Bernese which is just one of the GNSS data processing software. Bernese could generate regional ionosphere model as IONosphere map EXchange (IONEX) format by using Receiver Independent Exchange (RINEX) observation files cleaned during the Precise Point Positioning (PPP) process, and the estimation of ionosphere is made by setting the main parameter estimation program of Bernese (Dach et al. 2015). It is possible for Bernese to generate ionosphere by processing a RINEX observation file that was generated from a receiver which is not included in the CODE GIM. STEC is obtained by multiplying slant factor to this VTEC. Ionospheric delay in time domain is obtained by using conversion formula according to the signal frequency. Finally, ionospheric delay of the TWCP is calculated from Bernese and GIM.

The rest of this paper is organized as follows. The method of obtaining VTEC by GIM and Bernese is provided in Section 2. In Section 3, VTEC values generated by using Bernese are compared to those by GIM. In addition to this, ionospheric delay of the TWCP for Eutelsat 172A will be showed by considering ionosphere thin shell model with KRISS and Koganei (KGNI) stations. Finally, in Section 4, summarized conclusions are presented.

2. METHOD

Ionosphere reaches from about 50 km to 1000 km with vertical profile of the electron density according to the altitude, and the ionization is mainly occurred due to the extreme ultraviolet and X-ray emission (Komjathy 1997). The signals of the GPS satellites are refracted by electrons in the ionosphere, and this effect causes signal delays which are proportional to total electron content along the signal path (Datta-Barua 2008). Thin shell model, which assumes nearly all electrons are concentrated at a chosen altitude, is used for mapping TEC. TEC is represented in TEC units (TECU), and Eq. (1) explains how to convert TEC into ionospheric delay in time domain (Fujieda et al. 2014).

\(I=\frac{40.3\cdot TEC}{c\cdot f^2}\)                                                                                                                            

\(\textrm{where 1 TECU}=\ {10}^{16}\ \textrm{electrons}/m^2\)                                                                                  (1)

where \(I\) is ionospheric delay in time domain, \(c\) is the speed of light, is the signal frequency. For Eutelsat 172A, the uplink frequency is 14.314625 GHz, and the downlink frequency is 12.566625 GHz. In TWCP, ionospheric delay is included in the time difference formula between stations A and B as shown in Eq. (2) (Fujieda et al. 2014). In this analysis, A is the KRISS station, and B is the KGNI station which is IGS station near NICT.

\(\tau_a\left(t\right)-\tau_b\left(t\right) = {[ω+αt-ω-βt]+(2ω_u ω_d)I \over ω_+^2-ω_-^2}\)                                                                                          

\(\textrm{where} \ \ \ \omega_+\equiv\omega_u+\omega_d,\ \omega_-\equiv\omega_u-\omega_d\)                                                                                           

\(\alpha\left(t\right)\equiv\phi_{ab}\left(t\right)-\phi_{ba}\left(t\right),\ \beta\left(t\right)\equiv\phi_{aa}\left(t\right)-\phi_{bb}\left(t\right)\)                                                                               

\(I=[\left(I_{da}\left(t\right)-I_{ua}\left(t\right)\right)-\left(I_{db}\left(t\right)-I_{ub}\left(t\right)\right)]\)                                                                                  (2)

where \(\tau\left(t\right)\) is the times of the reference clocks at station,\(\omega_u\) and \(\omega_d\) is uplink and downlink signal angular frequency, \(\phi_{ij}\) is phase information from station i to j, and I is the ionospheric delay.

Fig. 1. Schematic diagram of spatial interpolation.

CODE GIM have latitude, longitude, and 1 hour resolution since 19 October 2014 (Dach et al. 2015). Spatial and temporal interpolation is applied to GIM in order to get VTEC at desired latitude, longitude, and time. The IONEX format and interpolation concept is explained by Schaer et al. (1998). Fig. 1 shows schematic diagram of the spatial interpolation, and Eq. (3) explains how to calculate. Nearest four VTEC are collected first, and then weighting is applied to the VTEC according to the distances. As the distance is more nearer, the more weight is counted.

\(E_s=\frac{1}{3}\frac{\left(R_s-R_1\right)}{R_s}E_1+\frac{1}{3}\frac{\left(R_s-R_2\right)}{R_s}E_2+\frac{1}{3}\frac{\left(R_s-R_3\right)}{R_s}E_3+\frac{1}{3}\frac{\left(R_s-R_4\right)}{R_s}E_4\)                                                          (3)

\(\textrm{where} \ \ R_s=\sum\limits_{i=1}^{4}R_i\)                                                                                                                             

where \(E_s\) is the spatial interpolated VTEC, \(E\) is the VTEC, and \(R\) is the distance between desired and each selected VTEC. Fig. 2 and Eq. (4) describes temporal interpolation method. After the spatial interpolation is performed, this VTEC is used for temporal interpolation. Since GIM has 1 hour time interval among VTEC maps, temporal interpolation is conducted on hourly maps by weighting more to the near time VTEC.

Fig. 2. Schematic diagram of temporal interpolation.

\(E_t\left(\beta,\lambda,t\right)=\frac{T_{i+1}-t}{T_{i+1}-T_i}E_{s(i)}\left(\beta,\lambda\right)+\frac{{t-T}_i}{T_{i+1}-T_i}E_{s(i+1)}\left(\beta,\lambda\right)\)                                                             (4)

\(\textrm{where} \ T_i \le t < T_{i+1}\)                                                                                                                 

where \(E_t\) is the temporal interpolated VTEC, \(E_s\) is the spatial interpolated VTEC, \(\beta\) is the latitude, is the longitude, and t is the time. Since orbit position almost stationary due to general GEO characteristics, one Eutelsat 172A position is used at 29 January 2017. This is obtained from Eutelsat 172A contact. IPP is defined as the location that the signal path reaches the altitude of the thin shell model (Datta-Barua 2008). Fig. 3 displays schematic diagram of the ionosphere thin shell model, and slant factor is applied from Eq. (5) (Schaer 1997). Since both GIM and Bernese uses 450 km as the height of the thin shell for the ionosphere, this value is also applied to the calculation of the IPP between station and Eutelsat 172A. Station VTEC is obtained from at the latitude and the longitude of the receiver. However, both GEO VTEC and GEO STEC are calculated from at the latitude and the longitude of the GEO IPP which is calculated from the geometry between receiver and GEO. Slant factor at IPP of the GEO is used to change GEO VTEC into GEO STEC.

Fig. 3. Ionosphere thin shell model.

\(GEO\ STEC=\ \frac{GEO\ VTEC}{\cos{z^\prime}}\)                                                                                                     

\(\textrm{with} \ \ z^\prime=\arcsin{\frac{R\cos{el}}{R+H}}\)                                                                                                      (5)

where \(z^\prime\) is the zenith angle at IPP, \(el\) is the elevation angle of the station, \(R\) is the mean radius of the earth, and \(H\) is the height of the thin shell model. From the geometry, earth-centered angle could be calculated from Eq. (6), and this is used to obtain the latitude and the longitude at the IPP.

\(\alpha=\frac{\pi}{2}-z^\prime-el\)                                                                                                               (6)

where \(\alpha\) is the earth-centered angle between station and a satellite, \(z^\prime\) is the zenith angle at IPP, and is the elevation angle of the station. Eq. (7) explains how to calculate the latitude and the longitude at the IPP (Datta-Barua 2008).

\(\beta_{IPP}=\arcsin{\left(n\right)},\ \lambda_{IPP}=\ \lambda_s+\arctan\left(\frac{e}{u}\right)\)                                                                                         

\(\textrm{where} \ \left(\begin{matrix}e\\n\\u\\\end{matrix}\right)=\ CBA\left(\begin{matrix}0\\0\\1\\\end{matrix}\right)\)                                                                                                             

\(A=\left(\begin{matrix}1&0&0\\0&\cos\alpha&\sin\alpha\\0&-\sin\ \alpha&\cos\ \alpha\\\end{matrix}\right),~~~ B = \left(\begin{matrix}\cos a& \sin a& 0\\-\sin a& \cos a& 0\\0&0&1\\\end{matrix}\right), ~~~ C=\left(\begin{matrix}1&0&0\\0&\cos{\beta_s}&\sin{\beta_s}\\0&-\sin{\beta_s}&\cos{\beta_s}\\\end{matrix}\right)\)                                                        (7)

where \(\beta_{IPP}\) is the latitude at the IPP, \(\lambda_{IPP}\) is the longitude at the IPP, \(\beta_s\) is the latitude of the station, \(\lambda_s\) is the longitude of the station, \(\alpha\) is the earth-centered angle, and is the azimuth angle of the station.

Bernese software version 5.2 is used to generate VTEC at desired latitude and longitude by using a station. Dach et al. (2015) explains ionosphere modeling and estimation at chapter 13. PPP DEMO of the Bernese Processing Engine is used to produce VTEC. The height of single layer is selected as 450 km, and the IONEX control file is changed to generate VTEC at desired latitude and longitude with 5 min time interval.

3. RESULTS

VTEC is compared at KRISS station location between CODE GIM and Bernese. Fig. 4 shows these VTEC with 5 min time interval at 29 January 2017. Two graphs show similar trends in this day, and this means that Bernese has generated VTEC appropriately. Bernese shows different variations between hours compared to the GIM around 19 UT. Since the ionosphere variations near 19 UT also showed in the raw data of the receiver from investigation, it is thought that the Bernese reflected those variations. In addition to this, the characteristics of the harmonic function, which is the ionosphere model of the Bernese, may be added to the ionosphere variations according to the setting process of the Bernese. According to the investigation on the representative influencing factors such as geomagnetic activity, earthquakes, and volcanoes, it is considered that the effects of geomagnetic activity and volcanoes seems to be small on the variations near 19 UT. It is necessary to study further on the influences of the earthquakes around the Korean Peninsula and other factors for accurate identification. VTEC difference between Bernese and GIM is shown in Fig. 5. The maximum difference is about 2.38 TECU in absolute value, and the Root Mean Square (RMS) is about 0.92 TECU. Since Bernese generate VTEC over the KRISS receiver using the RINEX observation file which is generated by the KRISS receiver, Bernese may be more suitable to represent the ionosphere over the KRISS than the GIM. In other words, CODE GIM could not accurately represent the ionosphere over the receiver of the KRISS because it only uses RINEX files from selected receivers around the world that does not include the receiver of the KRISS. The ionosphere over the receiver that used to generate the CODE GIM would be more accurate than that of the unused receiver. Although GIM gives relatively accurate VTEC, this could not overcome resolution in latitude, longitude, and time. On the other hand, Bernese could generate high resolution VTEC with 5 min time interval and first prime number in latitude and longitude by setting the IONEX control file. In this way, more precise IONEX file could be obtained by setting the IONEX control file in Bernese software compared to the IONEX of the GIM. In this perspective, Bernese would be suitable to generate high resolution ionosphere data. These gives confidence that Bernese is appropriate software to generate VTEC for regional analysis. Ya'acob et al. (2010) showed that the availability of the Bernese 5.0 in the generation of the TEC map over single receiver station by analyzing with CODE GIM. Abdullah et al. (2008) stated that Bernese 5.0 can be used to determine TEC over single station through the IGS station analysis with IGS GIM. Todorova et al. (2003) demonstrated that the regional ionosphere models from the Bernese 4.2 shows a good agreement with different sources such as GIM, very long baseline interferometry, and TOPEX/Poseidon satellite.

Fig. 4. VTEC over the KRISS station from Bernese and GIM on 29 January 2017.

Fig. 5. VTEC difference over the KRISS station between Bernese and GIM on 29 January 2017.

Position information of the KRISS and KGNI station are given in Table 1. Although the IPP is located somewhat distant from each stations, it is stated that single station could observe to the 1270 km on the condition that elevation cutoff angle is and single layer height is 450 km (Alcay 2016). Alcay (2016) has called this as coverage circle, and analyzed ionosphere variation related with earthquake by generating regional ionosphere map from single station with Bernese 5.0. These gives prove of adequacy to the Bernese software in the calculation of the VTEC at IPP.

Table 1. Position of the KRISS and KGNI.

Position	KRISS		KGNI
	Latitude (degree)	Longitude (degree)	Latitude (degree)	Longitude (degree)
Station IPP	36.4 32.8	127.4 134.1	35.7 32.3	139.5 143.8

 

Fig. 4 shows VTEC at the latitude and longitude of the KRISS receiver by processing KRISS data. In contrast, Fig.6 shows VTEC at the latitude and longitude of the IPP between KRISS and GEO by processing KRISS data. Similarly, Fig.7 shows VTEC at the latitude and longitude of the IPP between KGNI and GEO by processing KGNI data. The longitude difference between KRISS receiver and GEO IPP is , and this is corresponds to about 27 minutes in time theoretically. This difference was represented as the shift on the time axis, as shown in between Figs. 4 and 6. Figs. 6 and 7 shows VTEC at each station IPP, and they have similar variation trends between Bernese and GIM. In case KRISS, two methods show almost same VTEC, but they show significant differences at time interval around 19 UT as shown in Fig. 4. KGNI shows overall agreement in variation pattern. Fig. 8 shows VTEC difference between Bernese and GIM at each station IPP. KRISS has showed about 4.14 TECU in absolute value, and this is larger than about 2.62 TECU of the KGNI. This is mainly due to the time interval around 19 UT at KRISS. However, the RMS of the KRISS is about 1.03 TECU, and this is smaller than about 1.58 TECU of the KGNI. Fig. 9 shows VTEC difference between KRISS and KGNI at each IPP from Bernese and GIM. Bernese has showed maximum 3.3 TECU in absolute value with about 1.33 TECU at RMS, and these are larger than about 2.3 and 0.92 TECU of the GIM each. This difference between KRISS and KGNI would be represented as ionospheric delay in TWCP.

Fig. 6. IPP VTEC of the KRISS from Bernese and GIM on 29 January 2017.

Fig. 7. IPP VTEC of the KGNI from Bernese and GIM on 29 January 2017.

Fig. 8. IPP VTEC difference between Bernese and GIM from KRISS and KGNI on 29 January 2017.

Fig. 9. IPP VTEC difference between KRISS and KGNI from Bernese and GIM on 29 January 2017.

Ionospheric delay in TWCP is shown in Fig. 10. Bernese has showed from -4.77 to 5.69 ps, but GIM has showed from -1.53 to 2.54 ps. This means that ionospheric delay could reach up to 5.69 and 2.54 ps in Bernese and GIM each, and this indicates that results of the TWCP could be corrected by this amount. Fig. 11 shows ionospheric delay difference between Bernese and GIM in TWCP. The maximum is about 6.29 ps in absolute value, and the RMS is about 2.67 ps. This indicates that Bernese could correct the results of the TWCP with 6.29 ps rather than using GIM. Since VTEC at each station IPP was similar, the ionospheric delay in TWCP was also small. This may due to the short baseline distance between KRISS and KGNI. If the two region shows large different ionosphere states, ionosphere effect in TWCP would be large. Fujieda et al. (2014) stated that the measurement precision of the TWCP expected to be three orders of magnitude better than about 0.5 ns of the TWSTFT. Because of this high precision of the TWCP, the effect of the ionosphere has become significant. Therefore, TWCP has to remove the ionospheric delay rather than ignore that as in conventional TWSTFT. As for the TWCP link between NICT and Physikalisch-Technische Bundesanstalt, the ionospheric delay is over 100 ps (Fujieda et al. 2014). They explained that this is because of the local time difference and low elevation angle due to a very long baseline. In our case, the ionospheric delay of the Bernese for TWCP could be about 10 ps in peak-to-peak value. This is relatively small value rather than aforementioned, but it is still important to remove this amount of the ionospheric delay due to high precision of the TWCP. Furthermore, since the ionosphere of the KGNI could be changed significantly due to environmental effect like earthquakes, the ionospheric delay of the TWCP could be large.

Fig. 10. Ionospheric delay of the TWCP from Bernese and GIM on 29 January 2017.

Fig. 11. Ionospheric delay difference of the TWCP between Bernese and GIM on 29 January 2017.

4. CONCLUSIONS

In this paper, we presented a method of generating ionospheric delay for time comparison for a GEO satellite by using Bernese. The results has showed software adequacy with RMS as 0.92 TECU when compared with CODE GIM, and this has indicated suitability in the representation of the ionosphere for a specific receiver. Bernese has definite advantages in latitude, longitude, and time resolution compared to the GIM. At each station IPP, VTEC difference between KRISS and KGNI showed absolute maximum 3.3 and 2.3 TECU in Bernese and GIM, respectively. The overall VTEC differences of the two regions are small due to short baseline, and these amount are represented as ionospheric delay in TWCP. Ionospheric delay reached up to maximum 5.69 ps in Bernese, and the GIM has showed 2.54 ps. This means that Bernese could correct TWCP results with 5.69 ps. Moreover, Bernese could correct ionospheric delay up to 6.29 ps rather than using GIM in this day. Bernese is more suitable than CODE GIM in two perspectives. First, it generate ionospheric delay by processing RINEX observation file which is generated by the installed receiver at that location. Second, high resolution IONEX file could be obtained by setting the IONEX control file in Bernese compared to the IONEX file of the CODE GIM. However, CODE GIM also could be used for TWCP, and future studies could apply this approach using other software such as GIPSY-OASIS and GPSTk.

This paper has provided detailed description of how to obtain ionospheric delay between two locations, and this method is applied in a practical case. In this case study, the ionospheric delay of the Bernese is about 10 ps in peak-to-peak value. In atomic clock, the level uncertainty corresponds to approximately 10 ps for 1 day averaging time, therefore the time synchronization performance requires less than 10 ps. Currently, a satellite and ground station is synchronized within approximately 2ns, but the smaller the required performance, such as 1 ns or less, the better. Since the ionospheric delay could exceed over 100 ps in a long baseline unlike 10 ps in this relatively short baseline, it is considered important to compensate the ionospheric delay for more precise time synchronization of a satellite and ground control station. Although one specific case has been shown in this research, it has been shown that anyone could eliminate the ionospheric delay for the TWCP between two organizations at the receiver installed location. By using this technique, an organization could get high precision time comparison capability. Further, this high precision TWCP method could be used to compare atomic clocks, and also applied in the ground stations for the development of the future domestic satellite navigation system.