Derivation of TEC and estimation of instrumental biases from GEONET in Japan

This paper presents a method to derive the ionospheric total electron content (TEC) and to estimate the biases of GPS satellites and dual frequency receivers using the GPS Earth Observation Network (GEONET) in Japan. Based on the consideration that the TEC is uniform in a small area, the method divides the ionosphere over Japan into 32 meshes. The size of each mesh is 2 ◦ by 2 in latitude and longitude, respectively. By assuming that the TEC is identical at any point within a given mesh and the biases do not vary within a day, the method arranges unknown TECs and biases with dual GPS data from about 209 receivers in a day unit into a set of equations. Then the TECs and the biases of satellites and receivers were determined by using the leastsquares fitting technique. The performance of the method is examined by applying it to geomagnetically quiet days in various seasons, and then comparing the GPS-derived TEC with ionospheric critical frequencies ( foF2). It is found that the biases of GPS satellites and most receivers are very stable. The diurnal and seasonal variation in TEC and foF2 shows a high degree of conformity. The method using a highly dense receiver network like GEONET is not always applicable in other areas. Thus, the paper also proposes a simpler and faster method to estimate a single receiver’s bias by using the satellite biases determined from GEONET. The accuracy of the simple method is examined by comparing the receiver biases determined by the two methods. Larger deviation from GEONET derived bias tends to be found in the receivers at lower ( <30 N) latitudes due to the effects of equatorial anomaly.


Introduction
The total electron content (TEC) is one of the most important parameters used in the study of the ionospheric properties. Acting as a dispersive medium to the Global Positioning System (GPS) satellite signals, the ionosphere causes a group delay and a phase advance to the radio waves propagating from a GPS satellite to a ground-based receiver. TEC can be obtained from the difference in the group delays of dual-frequency GPS observations. However, there exists an instrumental delay bias in each signal of the two GPS frequencies. Their difference, referred to as instrumental or differential instrumental bias, affects the accuracy of the TEC estimation greatly. The combined satellite and receiver biases can even lead to a negative TEC.
The task of assessing GPS satellite and receiver biases has been assumption dependent and time consuming. Assuming that (1) the electron distribution lies in a thin shell at a fixed height above the Earth; (2) the TEC is time-dependent in a reference frame fixed with respect to the Earth-Sun axis; (3) the satellite and receiver biases are constant over several hours. Several authors (Lanyi and Roth, 1988;Coco et al., 1991) made their analysis with data from a single station during local nighttime, and they modeled the vertical TEC by a quadratic function of latitude and longitude. Wilson et al. (1992Wilson et al. ( , 1995 extended the thin spherical shell fitting technique to data sets from a GPS network in a 1-day or 12-h unit, and represented the vertical TEC as a spherical (surface) harmonic expansion in latitude and longitude. Sardòn et al. (1994) modeled the vertical TEC as a second-order polynomial in a geocentric reference system, where the coefficients of the polynomial are simulated with random walk stochastic processes. The coefficients (and hence, the TEC) and instrumental biases are then estimated by using a Kalman filtering approach. A common feature of the previous works is that an assumption of a rather smooth ionospheric behavior had to be introduced in the studies. Recently, with data collected from more than 1000 receivers of the GPS Earth Observation Network (GEONET) in Japan, Otsuka et al. (2002) produced two-dimensional maps of the TEC having a high spatial resolution of 0.15 • by 0.15 • in latitude and longitude. Although they removed the instrumental biases in order to derive the absolute vertical TEC, they did not discriminate between the satellite and receiver biases separately.
In this paper, we present a method to derive the TEC over Japan, and estimate the biases of GPS satellites and the dual P-code receivers that are part of GEONET in Japan. Our method is different from that of Otsuka et al. (2002) in that along with the TEC, both the satellite and the receiver biases can be obtained. The algorithm is depicted in detail in Sect. 2. We show in Sect. 3 the results of an application of the proposed method to three geomagnetically quiet days in the summer, autumn and winter of 2001, respectively. After the stability of the satellite biases is shown, day-to-day variation in instrumental bias is discussed. Evaluation of the GPSderived TEC is made by comparison with ionosonde's ionospheric critical frequency (referred to as foF2) observations. Discussion on the accuracy of the GEONET-based method is presented with the goodness of fit to the data. We propose in Sect. 4 a simpler and faster method to estimate a single receiver's bias by using its GPS observations and known satellite biases. The accuracy of the method is manifested by applying it also to the 9 days and by comparing the results with those in Sect. 3. The main results obtained are summarized in Sect. 5. Finally, the conclusions drawn are presented in Sect. 6.

TEC extraction from GPS observation
There are 28 GPS satellites currently orbiting the Earth at an inclination of 55 • and at a height of 20 200 km. They broadcast information on two frequency carrier signals, which are 15 7542 GHz (referred to as f 1 ) and 12 276 GHz (referred to as f 2 ), respectively. GPS observations give two distances (known as pseudorange) and two phase measurements corresponding to the two signals. Because of the dispersive nature of the ionosphere, the two radio signals are delayed by different amounts (known as group delay), and their phases are advanced when they propagate from a satellite to a receiver on the Earth. The slant path TEC sl from a satellite to a receiver can be obtained from the difference between the pseudoranges (P 1 and P 2 ), and the difference between the phases (L 1 and L 2 ) of the two signals (Blewitt, 1990) where k, related to the ionosphere refraction, is 80.62 (m 3 /s 2 ). λ 1 and λ 2 are the wavelengths corresponding to f 1 and f 2 , respectively. Because of the 2π Fig. 1. Geometry of a GPS satellite (S), the ionosphere, and a receiver (R). While the total electron content is retained, the ionosphere is assumed to be a screen sphere lying at the height of 400 km from the ground. Here, P represents the intersection of the line of sight and the ionosphere, χ is zenith angle.
ambiguity in the phase measurement, TEC sll from the differential phase is a relative value, but it has higher precision than TEC slp . To retain phase path accuracy for the slant path TEC sl , TEC sll are fitted to TEC slp , introducing a baseline, B rs , for the differential phase related TEC sll Horvath and Essex, 2000) TEC sl = TEC sll + B rs .
If having N measurements, the baseline B rs in this paper is computed as the average difference between pseudorangederived TEC slp i and phase-derived TEC slli over the index i from i = 1 to i = N inclusive.
where the square sine of the satellite's elevation α i is included as a weighting factor, as the pseudorange with low elevation angle is apt to be affected by the multipath effect and the reliability decreases. Consequently, the contribution to the baseline determination is greatly depleted from slant paths with low elevations. When making the above calculation of B rs , a data-processing step is included to identify possible cycle-slips in either L 1 or L 2 phase measurements (Blewitt, 1990). Thus, this study works with pseudorange-leveled carrier phases that are free of ambiguities and have lower noise and multipath effects than the pseudoranges. With a 30-s time series of dual GPS data, this part of the process is done for each pair of satellite receivers independently. All effects on the phases and pseudoranges that are common to both frequencies (such as distance of receiver satellite, clock offsets, tropospheric delay, etc.) of the obtained slant path TEC sl are removed, but frequency-dependent effects, like multipath and the differential instrumental biases in the satellite and the receiver, are still present. To convert to a vertical TEC from a slant path TEC sl , the ionosphere is assumed to be a thin screen shell encircling the Earth and its center is assumed to be the same as that of the Earth. The geometry of the GPS satellite, receiver and the ionosphere is shown in Fig. 1. The intersection of the slant path from the satellite (S) to the receiver (R) through the ionosphere is referred to as a piercing point (P ). The zenith angle χ is expressed as the following where α is the elevation angle of the satellite, R E is the mean radius of the Earth, and h is the height of the ionospheric layer, which is assumed to be 400 km in this paper. Further, setting satellite and receiver biases as b s and b r , respectively, then the vertical TEC is The determination of the absolute TEC and the instrumental biases will be described following an introduction of GEONET, a dense GPS receiver network in Japan.

GEONET in Japan and mesh division
GEONET is a GPS Earth Observation Network set up by the Geographical Survey Institute (GSI) of Japan. It has more than 1000 GPS receivers spread over Japan (Miyazaki et al., 1997), about 209 of which give precise code pseudoranges at both frequencies. As shown in Fig. 2, the nationwide distributed receivers form a sufficiently dense network, covering an area from 27 • N to 45 • N and from 127 • E to 145 • E in geographical latitude and longitude, respectively.
Also shown in the map of Fig. 2 are 32 meshes drawn with dashed lines, in which TEC should be evaluated independently. Each mesh is 2 • by 2 • in longitude and latitude, respectively. There are as many as 20 receivers in some of the meshes. There are several meshes with no receivers within. The TEC at these meshes can be obtained as well, because there are receivers in their adjacent meshes, and the piercing points spread widely depending on the satellite location and the numbers of satellites.

Determination of TEC and instrumental biases
Without employing a complex mathematical model, it is assumed in this study that the vertical is identical at any point within a mesh, but TECs for different meshes can differ. This means that the TEC is taken to be local time-independent within 8 min, if converting the mesh width of 2 • in longitude to local time. Hence, for those lines of sight converging on the same mesh, the vertical components of their slant path TECs are all taken to be the same. It is also assumed that the satellite and receiver biases do not vary within one day.
For the line of sight from satellite j to receiver k piercing through the ionosphere in mesh m at time t, referring to Eq. (6), we can write the following equation where i denotes the order of the measurement at time t. The unknowns in Eq. (7) are TEC i , b s j , and b rk . With 28 satellites, 209 receivers, using observations with 15 min interval, the absolute TEC at 32 meshes for one day, 3300 unknowns in total, can be estimated by solving the following set of equations expressed in matrices where the vector on the right-hand side consists of the slant path TEC sl . The number of the TEC sl in the vector is L. The vector on the left-hand side denotes unknowns of the TEC i , the satellite bias b sj , and the receiver bias b rk . The number of the unknowns is I + J + K. The matrix on the left-hand side of Eq. (8) consists of coefficients, secχ for TEC, 1 for b s , 1 for b r and 0. It has (I + J + K) × L elements. For one day, for each mesh there are 96 values of TEC, for 32 meshes the number of unknown TECs is 96 × 32, that is I = 3072; J = 28, representing 28 satellites; k = 209, being the receiver number. Because it is not possible to determine unambiguously all the satellites and receiver biases absolutely, one of them (normally one receiver) is set to be 0, as a reference.
Then with a least-squares fitting technique, the solution to the above set of equations can be obtained by the singular value decomposition (SVD), which avoids unrealistic solutions of the equation system (Press et al., 1992). In our practical calculation, the number of equations is about 35 000. It takes about 8 h to carry out the whole process, from reading the GPS data to solving the Eq. (8), by a personal computer (PC) using a Pentium 4 processor.

Results of an application of the method
In order to demonstrate the performance of the technique, several days around the solstice and equinox period of 15-17 June, 20-22 September, and 21-23 December 2001 were selected, before and during which it is geomagnetically quiet (K p < 4). With the procedure described above, instrumental biases and vertical absolute TEC over Japan for each day are obtained. The selected reference receiver is located at 34.16 • N, 135.22 • E, which has more than 10 receivers surrounding it in the same mesh. Figure 3 shows the estimated satellite biases for the 9 days over a six-month time span, as a function of the day of year. The vertical dashed lines divide the inconsecutive days.

Instrumental biases
Here, the biases are those relative to their means that are indicated in the lower part of the panel. For all the satellites each day, the mean of their biases is first computed, and this mean is then subtracted from each individual satellite bias (Coco et al., 1991).
Consequently, the systematic trends, such as changes in the reference receiver bias, have been removed from the satellite. Although the mean of the satellites biases decreased several ns (1ns = 2.853TECU, 1TECU = 2.853×10 16 e/m 2 ) from the summer to the winter, the relative biases are quite stable. Among satellite bias differences between inconsecutive days, even the largest value was about 1 ns. The standard deviation in bias was from 0.076 ns to 0.664 ns for the satellite biases for the 9 days. It is less than 0.5 ns for 19 of the 28 satellites. So, the day-to-day variation was very small for satellite biases.
The day-to-day variation of the estimated receiver biases was also small for most of the receivers. The distribution of the standard deviation of the receiver biases to the 9-day mean is shown in Fig. 4. The greatest value was about 4 ns. Sixty-nine percent of the receivers had a standard deviation in bias that was smaller than 1 ns; 93% had less than 2 ns. In Fig. 5, a scatter diagram relates the standard deviation in receiver bias for the 9 days to the geographical position of the receiver. It is evident that there is no latitude dependence of the receiver bias variation. This implies that ionospheric local characteristics have little effects on the instrumental bias determination. In spite of this, it is noticeable in Fig. 4 that there are several receivers (in mid-latitudes) with large dayto-day variation of biases. There might be several reasons for this, for example: (1) the unstableness in the receiver circuit itself; (2) bias variation of the reference receiver; (3) multipath effects. It is likely that the unstableness in the receiver is the most reasonable explaination, because the bias variation of the reference receiver would affect all the other receivers, and the multipath effects would not vary greatly day by day.  Fig. 4. Distribution of the standard deviation of the GEONET derived receiver biases from the 9-day mean. 93% of the cases are within 2 ns.

GPS-derived TEC
With the method described in Sect. 2, TEC over Japan can be determined at the same time as the instrumental biases. 15-min time series of TEC is shown in the top panel of Fig. 6 for the 9 days from the summer to winter of 2001, for a mesh at (35 • N, 139 • E). The vertical dashed lines separate inconsecutive days. In addition to diurnal features, seasonal variation is conspicuous. Data obtained by other observation techniques are useful for a verification of the GPS-derived TEC. Bottom-side sounding by ionosonde is operated routinely every 15 min at Kokubunji (35.7 • N, 139.5 • E). The value foF2, shown in the middle panel in Fig. 6, is used to evaluate the accuracy of the GPS-derived TEC. As is evident, the behavior of TEC is strikingly similar to that of the foF2. The variation in TEC and foF2 shows a high degree of conformity. This is also obvious for fine structures that are displayed in the daytime. These facts indicate that the GPS-derived TEC is mainly contributed from electrons in the F2-region. A more detailed comparison, the ratio of TEC to the square of foF2, is presented in the bottom panel of Fig. 6 for the 9 days. The diurnal and seasonal variation is clearly displayed. While the daytime level of the ratio is not much different from the summer to the autumn, it doubles in the winter, suggesting a greater contribution from the plasmaspheric electron content. Figure 7 shows contour maps of TEC over Japan in the  summer, the autumn and the winter of 2001. The TEC distribution has a simple pattern in the summer. The daytime TEC in the autumn has both a larger value and a larger gradient in latitude than that in the summer. It is even larger in the winter than that in the autumn. The nighttime TEC value in the winter is about half of that in the other two seasons.

Accuracy evaluation of the method
The standard deviation of the data from the fitting parameters (residuals) is used to measure how well the estimated parameters agree with the data (Bevington, 1969) where L is the number of the slant path TEC sl data (refer to Sect. 2.3). Table 1 lists the χ g values for the 9 days analyzed. χ g is less than 5 TECU for 7 days. It is about 8 TECU on 16 June 2001 (167). χ g is about 51 TECU on 22 September 2001 (265). Individual residual for each data point is examined for the day 265, on which χ g is extremely large. On this day the number of slant path TEC sl data used is 47 400. There are 12 991 data satisfying that   5, that is to say the fitting results agree well with most of the data. Furthermore, it is found that most of the large residuals are from those meshes at latitudes lower than 35 • ; among 1233 data yielding |TEC sljk − secχ jk TEC i − b sj − b rk |<10, there are 950 data are from meshes at latitudes lower than 35 • . It is probable that a steep latitude gradient in the low latitude ionosphere, created by the development of an equatorial anomaly in equinox, caused the large standard deviation in the fitting on day 265. Thus, the large residuals mainly come from the TEC gradient within meshes at lower latitudes. A large χ g , however, does not necessarily mean the low fitting accuracy of the instrumental biases; the estimated satellite biases on day 265 do not differ very much from those on day 264, as seen in Fig. 3. A comparison of the receiver biases on the two days is shown with a scatter plot in Fig. 8. The circles in the figure represent those receivers located at latitude 35 • , and the crosses refer to the receivers at latitudes ≤35 • . The agreement between the biases for the two days is very good, regardless of the receiver latitude, although moderate deviation can be found for a few receivers. Thus, even for the worst case in terms of residual, the method determines the instrumental biases with a high accuracy.

Estimation of bias for a single receiver
The method described in the above section is not always applicable to any situation, because the technique is based on a highly dense receiver network in a small area. Also, the algorithm requires a lengthy processing time, which does not meet the requirement of monitoring the ionosphere in nearly real-time. However, once the satellite biases are determined by using GEONET, those values can be commonly used in any other location in the globe, even where a single receiver is installed. This section will describe a simple and fast method to estimate the bias of a single receiver using the satellite biases determined by GEONET, and the accuracy of the simple method will be evaluated.

A simple method
Generally, one GPS receiver simultaneously receives signals from 5 or more GPS satellites at any time. The elevation angle of those satellites could vary widely. The piercing points would be scattered widely but within a limited area, roughly 23 • in longitude and 32 • in latitude, with the receiver at the center. From different satellites with different elevations the lines of sight to the receiver lead to a spatial variation of slant path TEC sl at any observation time. If the ionosphere is horizontally homogeneous and instrumental biases are correctly removed, the vertically converted TECs should be identical for all of the satellites. In an actual case, in which the ionosphere has a horizontal gradient and vertical structure, the scattering of vertical TECs is assumed to be the smallest when the instrumental biases are correctly removed. As the satellite biases are well determined by GEONET and shown to be stable (refer to Sect. 3), which are known values hereafter, the receiver bias is estimated independently from GEONET by trying out a series of bias candidates and finding the one that gives a minimum deviation of TECs to their mean. In a mathematical description, given a trial receiver bias b(i), the standard deviation of TECs to their mean is calculated at each observation time. Then, the total standard deviations, σ i , is obtained for the whole day. The value of b(i 0 ) when σ i takes the minimum value, σ i , is considered to be a correct receiver bias (hereafter, referred to as fitted receiver bias). It takes only several minutes to obtain the fitted receiver bias by a personal computer (PC) using a Pentium 4 processor.
When different receiver biases are applied, the dispersion of vertical TECs is examined by using actual data set. For the convenience of comparison, one receiver is chosen from GEONET, which is located at 35.53 • N, 137.89 • E. The results for the observations on 17 June 2001 are given in Fig. 9. The dashed lines are for slant path TEC sl from the satellites to the receiver. The solid lines represent vertically converted TECs after the satellite and receiver biases are removed. For the three panels, the satellite biases were identical and determined with the method described in Sect. 3, but the receiver bias was taken to be different: in the top panel, the receiver bias is a GEONET-derived one; in the lower two panels, the receiver biases were arbitrarily chosen so that it is much less than the GEONET-derived one in the middle panel, and much larger than the GEONET-derived one in the bottom panel. The corresponding value of σ i for each case is shown at the top right corner. It is evident that when an inappropriate  receiver bias is applied, the curves do not converge. Figure 10 shows the variation of σ i as a function of b(i) for the same data set. From the figure the receiver bias is determined as 2.78 ns, which is close to the value determined from GEONET, 2.29 ns. The difference between biases from the two methods is only 0.49 ns.

Accuracy of the simple method
The same procedure was applied to all the GEONET receivers, and the receiver biases derived from the two methods are compared. A scatter plot of the GEONET-derived bias versus the fitted bias on 17 June 2001 is shown in Fig. 11 for all receivers. The agreement between the GEONET b r and the fitted one is amazingly good. Figure 12 gives the distribution of the difference between the GEONET and the fitted biases, b r (= b r GEONET − b r f it ) (hereafter, refered to as an error of fitted bias or simply an error) for the same data set. It can be seen that for most of the receivers (93%), the errors are within ±2 ns. Table 2 summarizes the percentage of the number of receivers for which the errors are within ±2 ns for the 9 days analyzed. It is noticeable that on 22 September 2001 (the 265th day of the year) the fitted bias has a large error for about 1/3 of the receivers. Specifically, these receivers are located at latitudes lower than 35 • N, as shown in Fig. 13, where the error's latitude dependence for the other days is also displayed. This is in agreement with the large χ g on the day 265 discussed in Sect. 3.3. On the whole, the value of b r f it tends to be larger than that of b r GEONET for the receivers at lower latitudes (<30 • N), and the error tends to  increase with the decrease in latitude. This suggests that the ionospheric condition affects the bias determination by fitting for a single receiver. For further investigation of the error source, and hence, the limit in the application of the method, the total standard deviation of the TECs to their mean, σ , for each receiver was calculated by using the fitted receiver bias. The latitude variations of σ are shown in Fig. 14. By comparing Figs. 13 and 14, it can be seen that a large value of σ , or ill convergence, does not necessarily yield a large error. Taking 22 September 2001 as an example, the error decreased with the increase in σ at latitudes lower than 30 • N. The latitude dependence of the σ and hence, the bias error can be explained in terms of the TEC latitude gradient and the equatorial anomaly, which are clearly depicted in Fig. 14. Having high activity in the equinox, the equatorial anomaly is characterized by two electron density peaks (known as crest) in the vicinity of the geomagnetic latitude 15 • symmetric to the geomagnetic equator, which corresponds to about 25 • N geographically at Japan's longitude. For a receiver located at or near the crest of a equatorial anomaly, the satellites within the range tend to be distributed apart from the crest. The vertically converted TECs would have a mean smaller than the TEC through the crest. And the deviation of TECs from their mean, σ , would be smaller than that of TECs with large latitude gradient or variance.

Summary
The dual GPS data from 209 GEONET receivers in Japan was used to determine TEC over Japan, as well as the biases of satellites and receivers. The paper also proposed a faster and simpler way to estimate a single receiver's bias as long as the satellite biases are known. The methods described herein have been applied to geomagnetically quiet days in the summer, the autumn and the winter.
The main results obtained in the biases' estimation can be summarized as follows: 1. The standard deviation from the mean is from 0.076 ns to 0.664 ns for the 28 GPS satellite biases for 9 days over the six-month time span.
2. Ninety-three percent of the receiver biases have a standard deviation that is smaller than 2 ns from the mean for the 9 days. It can be as large as 4 ns for a few receivers.
3. The fitted bias for a single receiver is generally within Latitude dependence of difference between b r_GEONET and b r_fit Fig. 13. Latitude dependence of the difference between biases determined from the two different methods for the 9 days analyzed. The dashed line referring to no difference is plotted in each panel for easy comparison. from a GEONET derived bias tends to occur for those receivers at lower latitude (<35 • N) in the autumn and winter. This is the result from the steep latitude gradient in the local ionosphere, probably with the development of the equatorial anomaly effects.
Concerning the GPS-derived TEC, the following has been found from a comparison with foF2:

The diurnal and seasonal variations in TEC and foF2
show a high degree of conformity.
2. The ratio of TEC to the square of foF2 also showed diurnal and seasonal variation. The daytime peak value in the winter was about twice that in the summer and autumn.

Conclusions
It can be concluded based on the results of an analysis of data obtained from GEONET that the method described herein is efficient and qualified for use to derive the absolute TEC, and to determine the biases of GPS satellites and receivers. Since the day-to-day variation is small in satellite and receiver biases, it is only necessary that the instrumental biases be esti-mated or calibrated from time to time. This is especially true for satellite biases.
The proposed method for estimating a single receiver's bias is faster and sufficiently accurate for a receiver at midlatitude. It has the potential to meet the requirement of being able to monitor the ionosphere in nearly real-time. It can be also applied to the receiver far from a GPS network. But the accuracy of a fitting bias can be low for a receiver at a lower latitude, due to the effects of equatorial anomaly. This disadvantage can be avoided by determining the receiver bias at mid-latitude before its establishment at a lower latitude.
The GPS-derived TEC is mainly contributed from the electrons in the F2-region. It is shown from the ratio of TEC to the square of foF2 that plasmaspheric electron content is larger in the winter than that in the summer or autumn.