Estimation of GPS Differential Code Biases Based on Independent Reference Station and Recursive Filter

The differential code bias (DCB) of the Global Navigation Satellite Systems (GNSS) receiver should be precisely corrected when conducting ionospheric remote sensing and precise point positioning. The DCBs can usually be estimated by the ground GNSS network based on the parameterization of the global ionosphere together with the global ionospheric map (GIM). In order to reduce the spatial-temporal complexities, various algorithms based on GIM and local ionospheric modeling are conducted, but rely on station selection. In this paper, we present a recursive method to estimate the DCBs of Global Positioning System (GPS) satellites based on a recursive filter and independent reference station selection procedure. The satellite and receiver DCBs are estimated once per local day and aligned with the DCB product provided by the Center for Orbit Determination in Europe (CODE). From the statistical analysis with CODE DCB products, the results show that the accuracy of GPS satellite DCB estimates obtained by the recursive method can reach about 0.10 ns under solar quiet condition. The influence of stations with bad performances on DCB estimation can be reduced through the independent iterative reference selection. The accuracy of local ionospheric modeling based on recursive filter is less than 2 Total Electron Content Unit (TECU) in the monthly median sense. The performance of the recursive method is also evaluated under different solar conditions and the results show that the local ionospheric modeling is sensitive to solar conditions. Moreover, the recursive method has the potential to be implemented in the near real-time DCB estimation and GNSS data quality check.


Introduction
Nowadays, Global Navigation Satellite Systems (GNSS) observations provide various ways to estimate geophysical parameters, and one of the most important parameters is the total electron content (TEC) measurements for Earth's ionosphere research [1][2][3][4][5][6]. Due to the frequency dispersive property of the ionosphere, the ionospheric delay experienced by electromagnetic signals can be estimated by dual-frequency measurements. The relatively high-precision ionospheric TEC can be derived from dual frequency GNSS carrier phase leveling pseudorange measurements, in which the differential code bias (DCB) is one of the main errors that cannot be ignored for absolute TEC estimation. 2 of 17 Recently, a number of algorithms for GNSS DCB estimates based on multi-frequency measurements have been proposed [7][8][9]. Apart from the commonly adopted method by setting the DCB as a constant during the GNSS TEC estimation [10][11][12][13][14][15], some other algorithms are proposed to decrease the computation costs by utilizing global ionospheric maps [4]. Furthermore, various optimization methods have been proposed for faster DCB determination and more accurate ionospheric modeling, such as optimization on DCB estimation based on regional or single station ionospheric modeling [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Keshin [30] estimated the single receiver DCBs using the global ionospheric map (GIM) vertical TEC gridded values, based on least square techniques with a linear constraint. An algorithm called IGG algorithm based on generalized triangular series function and satellite filtering was developed and used to generate the DCB product from the Chinese Academy of Sciences [31]. Sarma et al. [32] proposed an algorithm based on singular value decomposition to estimate three receiver biases. Among all these methods, several fundamental assumptions are adopted. One of those assumptions is that all the electrons in the earth's ionosphere are concentrated in one single thin layer. While these single-layer TEC models may adopt different ionospheric effective heights, ranging from 350 to 450 km according to the spatial and temporal variation of the electron maximum height of the F2 layer. Another assumption is that both satellite and receiver DCBs are constant during a single day. One more constraint to overcome the ill-posed problem of the design matrix is usually called zero-mean constraint, which is used as a Lagrange multiplier to separate the satellite and receiver DCBs into two parts. It is necessary to point out that the zero-mean condition can vary with respect to the observable satellites, which may lead to a drift in DCB estimates with respect to DCB products provided by some IGS analysis centers like Center for Orbit Determination in Europe (CODE). In addition, Hong et al. [33] proposed a new and efficient algorithm using the geometry conditions between satellite and tracking receivers to estimate the receiver DCBs with no use of the single layer assumption, but requiring to know one receiver DCB.
For single-station DCB estimates, there are two basic approaches. One category of research models is based on a polynomial of coordinates in the Earth-Sun reference system, including the satellite and receiver DCBs in the model. The polynomial coefficients and DCBs are regarded as unknown parameters to be solved, and the observations form a linear set of equations which can be solved in the least square sense [34][35][36][37][38][39][40]. The other category makes use of the assumption that VTECs computed from different satellites over a certain angle of elevation are close to each other [41,42]. Therefore, the receiver DCBs can be estimated by minimizing the standard deviation of VTEC. Both of these two methods can be used to estimate the receiver DCBs for one single station, as long as the satellite DCB can be provided by other sources, such as the global GNSS network.
The ionospheric conditions vary with many factors such as space weather, local time and solar activities, geographic location and so on, particularly in equatorial regions in modeling and computation of TEC [43]. In order to estimate precise DCBs and analyze the performance of each tracking station located at different regions on the local ionospheric model, in this paper we present a recursive method for both satellite and receiver DCB estimates, together with a reference station selection procedure. The influence of the reference station selection on DCB estimates is discussed through a comparison between using all stations and only reference stations. The accuracy of DCB estimates and local ionospheric models are evaluated by the CODE DCB product and GIM. In the following, methods and observations are presented in Section 2, results and analyses are shown in Section 3, some discussions are presented in Section 4 and finally conclusions are given in Section 5.

Carrier Phase Smoothing Pseudorange
The time delays of GNSS signals received by GNSS receivers are converted to pseudo-range values and the phase shifts are recorded as phase delays in the receivers [44]. The standard model for pseudo-range and phase recordings for dual frequencies f 1 and f 2 are as follows [43]: where the subscript u denotes the receiver, the superscript m denotes the satellite, the subscript i denotes the frequency, p is the actual range between satellite and receiver, δt u and δt m are the clock errors for the receiver and satellite, respectively, d m trop,u denotes the troposphere group delay, d m ion i,u denotes the ionospheric delay, c is the speed of light in vacuum and DCB s denotes the satellite differential code bias (DCB) and DCB r denotes the receiver DCB, λ i is the wavelength, ϕ m trop,u denotes the phase shift due to the troposphere, ϕ m ion i,u denotes the phase shift due to the ionosphere and N m i denotes the initial phase ambiguity.
The slant TEC (STEC) can be obtained by utilizing the geometry-free observations with ignoring the higher order effects of the ionosphere. A relatively accurate TEC value can be obtained by a carrier phase smoothing pseudorange method, which can be expressed as follows: where L i is the carrier phase observation, N i is the carrier phase ambiguity, P i is the pseudorange observation, ε L and ε P are the noise and multi path errors in the carrier phase observation and pseudorange observation, respectively, f i is the frequency of L i , DCB r is the LEO receiver DCB, DCB s is the Global Positioning System (GPS) satellite DCB, and N is the average ambiguity.

Local Ionospheric Modeling
Instead of using the Global Ionospheric Map (GIM) as a reference ionosphere background or global ionospheric modeling, we establish a local ionospheric VTEC model centered at one single station based on a recursive filter to estimate the GPS satellite and ground receiver DCBs. The measurements of TEC plus DCB, which is denoted as relative STEC, can be expressed as follows.
where STEC rel represents the relative STEC, i.e., the combination of satellite and receiver DCB, STEC represents the absolute STEC, and ε represents all possible errors. Our local ionosphere model parameterizes the vertical TEC (VTEC) distribution by a single layer at height 450 km. The spatial-temporal variations of the VTEC near the ground GNSS receiver are parameterized by a 2-degree polynomial fitting.
where the five model coefficients are a 0 , a 1 , a 2 , a 3 , a 4 and a 5 . The coordinates longitude u 1 and latitude u 2 are both evaluated at the corresponding ionospheric pierce point (IPP) in a local sun-fixed coordinate Remote Sens. 2020, 12, 951 4 of 17 system, centered at the receiver position. The slant ionospheric delay is converted to the vertical ionospheric delay by an elevation-dependent mapping function MF.
where VTEC represents the vertical TEC at the ionospheric pierce point (IPP), MF represents the mapping function, θ represents the elevation angle, R e represents the mean radius of the Earth and H ion represents the single layer height, i.e., 450 km.

Model Initialization and Propagation
In order to maximize the reliability of the local VTEC model, we minimize the cost function at the initial state, with the observations occurring from 0000 LT to 0400 LT, due to the fact that the nighttime ionosphere is less variant than the daytime ionosphere.
The cost function is given by where x n represents the local ionospheric model parameters, d n represents the combination of satellite and receiver DCBs, called as combination for brevity later, A n represents a design matrix which is divide into two parts A nx and A nd , y represents the relative STEC, C 0 represents the initial covariance matrix, θ represents the elevation angle and σ relTEC represents the standard deviation of the relative STEC provided by the pre-processing procedure. Subscript "n" represents the parameter evaluated at step n when n is larger than 0. Note that the observable GNSS satellites between two adjacent steps are not necessarily the same. The elevation mask is set to 20 degrees to reduce the multipath effects and mapping errors. At the initial step, observations, of which the longitude differences between the corresponding IPP and model origin (0200 LT) are more than 30 degrees (2 h), are discarded.
After initializing the local ionospheric model, we proceed to propagate the model into the next step, assuming that the only difference of local ionospheric models between two adjacent periods results from the second term in (5), i.e., u 1 . The propagation procedure when neglecting the second order contribution can be expressed as follows: where x ' n represents the ionospheric model parameter vector before updating, x n−1 represents the ionospheric model parameters after updating at n − 1th step, U represents the model propagation matrix and δt represents the time interval between nth step and n − 1th step, which is set to 900 s in this study.
In terms of covariance matrix determination, it is assumed here that the model updating error is independent of the model error, which allows us to determine the total covariance in the following step by C n = C M n + C U n−1 (9) where C n represents the total covariance matrix, C M n represents the model covariance matrix and C U n−1 represents the model update covariance matrix. According to the model propagation procedure (8), the model covariance matrix can be expressed simply as where C M n represents the model covariance matrix at nth step and C M n−1 represents the model covariance matrix at n − 1th step. Furthermore, the model updating uncertainty is attributed mainly to the asymmetric term in the model propagation matrix U, i.e., a 1 δu 1 , due to the fact that the model propagated from the previous step may not be the optimal local ionospheric model in the current step. Thus, we assume that the model updating matrix is proportional to a 2 1n . The model updating covariance matrix can be written as follows: where a 1n represents the model coefficient a 1 in step n, C U n represents the model updating covariance matrix, C posterior n represents the posterior covariance matrix after updating the model, C ' n represents the posterior covariance matrix before updating the model, P represents the weight matrix determined by the inverse matrix of C 0 in (7),x n represents the model parameter vector after updating, N represents the number of all available observations, y n represents the relative STEC in step n and Y n represents the absolute STEC vector calculated from the ground observations. Note that in (12) and (13), we assume that all DCBs are constant or much less variant than ionospheric VTEC during every consecutive observation series, so that the variance evaluated in (12) is exactly the ionospheric model variance.

Recursive Filter
Considering the contribution of the local ionospheric model, the extended cost function can be expressed as To avoid any negative values of a 0 , which violates fundamental physics, we here employ an inequality constrained least square (ICLS) technique on the minimization process (15) [45]. Furthermore, the model propagation always continues in the same way, while we discard the DCB combinations derived from those observations, in which the first model parameter a 0 is smaller than 0 total electron content unit (TECU), which implies possibly poor qualities among them. After deriving the DCB combinations, equations containing satellite and receiver DCBs can be written in the form as (16).
where D is the design matrix, the elements corresponding to the satellite and receiver between which the observation happens are 1 and 0 for others, ε r is the residual, and Z is the combination vector of satellite and receiver DCBs estimated above. The satellite and receiver DCBs can be separated by imposing a zero-mean constraint on all the satellite DCBs, which can be expressed as: where S is the zero-mean constraint matrix, the elements corresponding to the satellite DCBs are 1 and 0 for others,X is the optimal solution for satellite and receiver DCBs, G is the weight matrix for the combination, σ com is the standard deviation vector for the combination evaluated separately in every consecutive observation series, N 0 is the total number of the DCB combinations and σ is the weighted residual. See the flowchart of the recursive method in Figure 1 for more illustrations.
Remote Sens. 2019, 12, x FOR PEER REVIEW 6 of 17 the combination, is the standard deviation vector for the combination evaluated separately in every consecutive observation series, is the total number of the DCB combinations and is the weighted residual. See the flowchart of the recursive method in Figure 1 for more illustrations.

Reference Station Selection
Many previous researches suggested that data provided by mid-latitude stations are more reliable and recommended to be used than those provided by low-latitude stations. Since the quality of data recorded by the IGS stations are not evaluated in pre-processing, it is necessary to discard some stations with worse performances compared with other stations in order to ensure the accuracy of ionospheric modeling, as well as the DCB estimates. An automatically iterative process is employed to select the reference stations, which can be divided into three steps as following: 1. Estimating the satellite and receiver DCBs based on the recursive method described above with all observations.

Reference Station Selection
Many previous researches suggested that data provided by mid-latitude stations are more reliable and recommended to be used than those provided by low-latitude stations. Since the quality of data recorded by the IGS stations are not evaluated in pre-processing, it is necessary to discard some stations Remote Sens. 2020, 12, 951 7 of 17 with worse performances compared with other stations in order to ensure the accuracy of ionospheric modeling, as well as the DCB estimates. An automatically iterative process is employed to select the reference stations, which can be divided into three steps as following: 1.
Estimating the satellite and receiver DCBs based on the recursive method described above with all observations. 2.
Estimating only the receiver DCBs with the initial satellite DCB estimates in the first step, by using the same recursive method. Note that the receiver DCBs in this step are determined by the median value of all the 15-min estimates, because in this step there is no ill-posed problem like in step 1.

3.
Comparing the receiver DCBs in the first and second step. If the difference between two receiver DCBs exceeds a threshold, we remove the corresponding station out of the set of reference stations and go back to step 1, until no station is removed in step 3. Then extract the final DCB estimates in this step. The threshold can be determined by the value of σ in the last formula in (17).

Experimental Data
In order to test the feasibility of the recursive filter method and reference station selection procedure, two experimental cases are designed. In each case, the GPS data are gathered from 130 IGS stations covering the whole globe (25 in low latitude region; 85 stations in mid latitude region; 20 stations in high latitude) from 15th to 31th January 2011. To avoid any extra uncertainties from a variable number of stations, we select those stations which provided continuous observations during the period of interests. In the first case, all stations play the equivalent role in in the DCB estimates. In the second case, however, those stations, which may lead to instabilities of DCB estimates or inaccuracies of ionospheric modeling, are discarded from the set of reference stations according to the criterion in step 3 mentioned above. The comparisons of DCB estimates between case 1 and 2 show the influence of the station selection on the performance of the recursive method. Moreover, the weight of each observation is evaluated as the second expression in (5). According to the daily DCB estimation of the recursive method described in Section 2, the combination of satellite and receiver DCBs is estimated in each time interval first and then the satellite and receiver DCBs are separated based on the zero-mean condition. It is noted that the discarded receiver's DCBs are considered as the results in the last estimation round, in which this receiver is not discarded by then. As is mentioned in Section 2.3, local ionospheric modeling is based on the observations occurring between 0000 LT and 0400 LT, in order to minimize the initialization error once per local day. Data in two consecutive days in universal time are needed and gathered together as a whole observation in one local day, station by station. However, the CODE receiver DCBs are estimated once per day in universal time rather than in local time. Thus, we have to align our receiver DCB results to CODE DCB by linear interpolation.
where DCB n is the aligned DCB estimate of day n in the universal time, lon is the longitude of the station in rad, DCB n−1 , DCB n and DCB n+1 are the original DCB estimate of day n − 1, n and n + 1 in the local time, respectively. Note that lon ranges from −π to π corresponding to 180 • W and 180 • E. It is also necessary to note that the satellite DCBs here do not need to be aligned due to the fact that satellite DCBs are very stable during consecutive days without hardware operations. Considering the shift between the universal time and local time in each station, the results in the last day, i.e., 31st January, are not shown below in the final estimates. Figure 2 shows the comparisons between the VTEC derived from the recursive method and interpolation using CODE GIM. The VTEC derived from the recursive method is generally consistent with that of CODE GIM. However, considering that the interpolation could result in considerable errors and the spatial-temporal resolution of GIM is worse than the that of single-station VTEC, there could possibly exist a discrepancy between GIM and recursive VTEC. Figure 3 shows the second and third parameters representing zonal (blue line) and meridional (red line) gradients of VTEC over GODZ station. We can see from Figure 3 that the sign of meridional gradients shows an obvious hemispheric difference. For the GODZ station, which is located at the northern hemisphere, the meridional gradients stay negative over the course of the day. This can be expected, since the VTEC over equatorial regions are higher than those over higher latitudinal regions. For zonal gradients, the results show diurnal variations, which means that the zonal gradient increases in the first half of the day (or more precisely before 1400 LT) and then decreases in the second half of the day, considering the ionization peak at around 1400 LT. Figure 4 shows the semi-monthly mean difference and average RMS, with respect to interpolated CODE GIM over 20

Reference Station Selection
As mentioned earlier, some stations which performed worse than the other stations during the periods of interest can be removed from the filter. Since more mid-latitude stations are removed from the set of reference stations compared with other days, we choose to show the abandoned stations on 17th January as the red dots in Figure 5. In order to analyze the spatial distribution of discarded stations, we divide one hemisphere into three parts, i.e., low (0-30 degrees), mid (30-60 degrees) and high latitude (60-90 degrees). Table 1 shows the variations of the latitudinal distributions of discarded stations during the period of interests. From Table 1 we can see that the discarded stations are mostly located at low and mid latitude region, which is expected, because the number of stations located at mid latitude is much larger than those located at low and high latitude and previous work has shown that the stations in low latitudes are generally not as good as those in mid latitudes [31]. However, it can be found that the number of discarded stations located at low latitude is much larger than the number of discarded stations located at mid latitude, even though there are only five more stations located in mid latitude. It can be concluded that stations located in the low latitude region generally performed worse than those located in the mid and high latitude region. Figure 6 shows the time series of GPS satellite DCBs estimated by the recursive method before (a) and after (b) selecting reference stations during the period of interests. In the red ellipsoid of Figure 6a, an apparent jump in satellite DCB of PRN 20 occurred in 17th January. After selecting reference stations, the negative influence from those bad stations can be mitigated according to the recovery of the jump. The reason is the sudden increase of the number of discarded stations in the mid latitude region on 17th January. In order to illustrate more explicitly the variation of distribution of discarded stations, we present the successive daily changes of geographic distribution of the reference stations (shown in green dot) and discarded stations (shown in red dot), during 15th to 30th in January 2011 in the supplementary materials Figure S1. As mentioned earlier, some stations which performed worse than the other stations during the periods of interest can be removed from the filter. Since more mid-latitude stations are removed from the set of reference stations compared with other days, we choose to show the abandoned stations on 17 th January as the red dots in Figure 5. In order to analyze the spatial distribution of discarded stations, we divide one hemisphere into three parts, i.e., low (0-30 degrees), mid (30-60 degrees) and high latitude (60-90 degrees). Table 1 shows the variations of the latitudinal distributions of discarded stations during the period of interests. From Table 1 we can see that the discarded stations are mostly located at low and mid latitude region, which is expected, because the number of stations located at mid latitude is much larger than those located at low and high latitude and previous work has shown that the stations in low latitudes are generally not as good as those in mid latitudes [31]. However, it can be found that the number of discarded stations located at low latitude is much larger than the number of discarded stations located at mid latitude, even though there are only five more stations located in mid latitude. It can be concluded that stations located in the low latitude region generally performed worse than those located in the mid and high latitude region. Figure 6 shows the time series of GPS satellite DCBs estimated by the recursive method before (a) and after (b) selecting reference stations during the period of interests. In the red ellipsoid of Figure 6a, an apparent jump in satellite DCB of PRN 20 occurred in 17 th January. After selecting reference stations, the negative influence from those bad stations can be mitigated according to the recovery of the jump. The reason is the sudden increase of the number of discarded stations in the mid latitude region on 17 th January. In order to illustrate more explicitly the variation of distribution of discarded stations, we present the successive daily changes of geographic distribution of the reference stations (shown in green dot) and discarded stations (shown in red dot), during 15th to 30th in January 2011 in the supplementary materials Figure S1.      Figure 7 shows the number of total numbers of DCB combinations by bar plot in case 1 (red) and in case 2 (green), while the weighted residual in Eq. (17) is presented by line plot in case 1 (blue) and case 2 (black). We can see that the total number of DCB combinations in case 1 stays very invariant in the period of interest, whilst that in case 2 varies from day to day. Comparing the red and green bars, decreases by 40% in the mean sense in case 2. We have to note though, that on 17 th January, is even larger than many other days; from Table 1 below, the number of discarded stations is largest during this period. This fact verifies that we correctly discarded those stations with only few and interrupted observations and kept the stably and continuously tracking stations in the set of reference stations. By comparing the blue and black lines, we can find that in most days except for 19 th and 31 st January, the weighted residual decreases in case 2, by a mean value of 0.10 ns. The comparison between two cases further confirms the feasibilities of the proposed reference selection procedure.  Figure 7 shows the number of total numbers of DCB combinations N 0 by bar plot in case 1 (red) and in case 2 (green), while the weighted residual σ in Eq. (17) is presented by line plot in case 1 (blue) and case 2 (black). We can see that the total number of DCB combinations in case 1 stays very invariant in the period of interest, whilst that in case 2 varies from day to day. Comparing the red and green bars, N 0 decreases by 40% in the mean sense in case 2. We have to note though, that on 17th January, N 0 is even larger than many other days; from Table 1 below, the number of discarded stations is largest during this period. This fact verifies that we correctly discarded those stations with only few and interrupted observations and kept the stably and continuously tracking stations in the set of reference stations. By comparing the blue and black lines, we can find that in most days except for 19th and 31st January, the weighted residual σ decreases in case 2, by a mean value of 0.10 ns. The comparison between two cases further confirms the feasibilities of the proposed reference selection procedure.
Remote Sens. 2019, 12, x FOR PEER REVIEW 11 of 17 Figure 7. The left axis represents the total numbers of DCB combinations (bar plot) in case 1 (red) and case 2 (green) and the right axis represents the weighted residuals (line plot) in Equation (17) in case 1 (blue) and case 2 (black). The horizontal axis represents the day of year in 2011.

Comparision with CODE DCB Product
In order to check the stabilities and accuracy of DCB estimates based on the recursive method, we introduce the day scatter and root mean square errors compared to CODE. The day scatter represents the day-to-day variation in daily DCB estimates, with respect to the monthly mean DCB estimates, as defined in (19) [45]. RMS is the root mean square errors with respect to CODE DCB product, which is expressed as (20).
where DS is the day scatter, D is the number of days of interests, is the satellite DCB estimates determined by recursive method, is the mean value of the satellite DCB estimates, is the satellite DCB from CODE DCB product. Figure 8 shows the day scatter of satellite DCB estimates determined by recursive method after selecting reference stations and corresponding RMS values with respect to CODE DCB during the periods of interest. The largest value of day scatter is 0.15 ns for PRN 19, while the lowest value is 0.05 ns for PRN 10; The largest value of RMS is 0.28 ns for PRN 17, while the lowest value is 0.05 ns for PRN 2. Figure 9 shows the difference of RMS and day scatter before and after selecting reference stations during the periods of interest. Comparing the results from two experimental cases, it can be seen that most of the satellite DCB estimates become more stable after selecting the reference stations, especially for PRN 20. The day scatter and RMS value for PRN 20 decreased by 0.10 ns and 0.11 ns, respectively. It can be suggested that the DCB stabilities of only the satellites that discarded stations may be able to track can be effectively improved. This is the reason why DCBs of satellites which cannot be tracked by the discarded stations do not change significantly after choosing the reference stations. For those satellites of which the day scatters increase, the increments of day scatters are all less than 0.03 ns, which is a slight increase compared with the decreases of some satellites, such as PRN 16, 20 and 27. Thus, we conclude that selecting reference stations are generally beneficial for the satellite DCB estimates.

Comparision with CODE DCB Product
In order to check the stabilities and accuracy of DCB estimates based on the recursive method, we introduce the day scatter and root mean square errors compared to CODE. The day scatter represents the day-to-day variation in daily DCB estimates, with respect to the monthly mean DCB estimates, as defined in (19) [45]. RMS is the root mean square errors with respect to CODE DCB product, which is expressed as (20).
where DS is the day scatter, D is the number of days of interests, B d is the satellite DCB estimates determined by recursive method, B is the mean value of the satellite DCB estimates, B d CODE is the satellite DCB from CODE DCB product. Figure 8 shows the day scatter of satellite DCB estimates determined by recursive method after selecting reference stations and corresponding RMS values with respect to CODE DCB during the periods of interest. The largest value of day scatter is 0.15 ns for PRN 19, while the lowest value is 0.05 ns for PRN 10; The largest value of RMS is 0.28 ns for PRN 17, while the lowest value is 0.05 ns for PRN 2. Figure 9 shows the difference of RMS and day scatter before and after selecting reference stations during the periods of interest. Comparing the results from two experimental cases, it can be seen that most of the satellite DCB estimates become more stable after selecting the reference stations, especially for PRN 20. The day scatter and RMS value for PRN 20 decreased by 0.10 ns and 0.11 ns, respectively. It can be suggested that the DCB stabilities of only the satellites that discarded stations may be able to track can be effectively improved. This is the reason why DCBs of satellites which cannot be tracked by the discarded stations do not change significantly after choosing the reference stations. For those satellites of which the day scatters increase, the increments of day scatters are all less than 0.03 ns, which is a slight increase compared with the decreases of some satellites, such as PRN 16, 20 and 27. Thus, we conclude that selecting reference stations are generally beneficial for the satellite DCB estimates.  For investigating the influence of reference station selection on the receiver DCB estimates, the comparisons of the differences between sample receiver DCBs with CODE receiver DCBs in the two cases mentioned above are presented in Table 2. From Table 2, we can see that most of the receiver DCBs, except for GODZ and QAQ1, become more stable in case 2 than in case 1, while the stabilities of receiver DCB determined by the recursive method are slightly worse than CODE DCB product, possibly because we use the observations with GPS satellites only. The monthly mean differences, with respect to CODE receiver DCB, decrease by 0.18, 0.23, 0.22, and 0.19 ns for HERS, IQAL, PERT and TIXI in case 2, while the other receiver DCBs remain at the same accuracy level. Considering the interactions among stations through the zero-mean condition, the reason may be attributed to that the discarded stations have little intersected satellites tracked by GODZ, QAQ1, REYK and YSSK. The influence of reference station selection on one station depends on the number of observations recorded, with the intersected satellites tracked by those discarded stations. Table 2. Monthly mean difference with respect to CODE receiver DCBs (MAD) and day scatter (DS) of eight sample stations during the second half of January 2011. Note that the number denotes  For investigating the influence of reference station selection on the receiver DCB estimates, the comparisons of the differences between sample receiver DCBs with CODE receiver DCBs in the two cases mentioned above are presented in Table 2. From Table 2, we can see that most of the receiver DCBs, except for GODZ and QAQ1, become more stable in case 2 than in case 1, while the stabilities of receiver DCB determined by the recursive method are slightly worse than CODE DCB product, possibly because we use the observations with GPS satellites only. The monthly mean differences, with respect to CODE receiver DCB, decrease by 0.18, 0.23, 0.22, and 0.19 ns for HERS, IQAL, PERT and TIXI in case 2, while the other receiver DCBs remain at the same accuracy level. Considering the interactions among stations through the zero-mean condition, the reason may be attributed to that the discarded stations have little intersected satellites tracked by GODZ, QAQ1, REYK and YSSK. The influence of reference station selection on one station depends on the number of observations recorded, with the intersected satellites tracked by those discarded stations.  For investigating the influence of reference station selection on the receiver DCB estimates, the comparisons of the differences between sample receiver DCBs with CODE receiver DCBs in the two cases mentioned above are presented in Table 2. From Table 2, we can see that most of the receiver DCBs, except for GODZ and QAQ1, become more stable in case 2 than in case 1, while the stabilities of receiver DCB determined by the recursive method are slightly worse than CODE DCB product, possibly because we use the observations with GPS satellites only. The monthly mean differences, with respect to CODE receiver DCB, decrease by 0.18, 0.23, 0.22, and 0.19 ns for HERS, IQAL, PERT and TIXI in case 2, while the other receiver DCBs remain at the same accuracy level. Considering the interactions among stations through the zero-mean condition, the reason may be attributed to that the discarded stations have little intersected satellites tracked by GODZ, QAQ1, REYK and YSSK. The influence of reference station selection on one station depends on the number of observations recorded, with the intersected satellites tracked by those discarded stations.

Dependence on Solar Condition
Considering that the ionosphere can be affected by solar activities severely, the performance of DCB estimates and ionospheric modeling can also degrade under disturbed solar conditions [46]. In order to evaluate the performance of recursive methods under different solar conditions, we carry out the same experiment described in Section 2 during 3rd to 12th in May 2015, when F10.7 index underwent an increase from 100 to 160 flux units, while F10.7 index kept a stable value of about 80 flux units during January 2011. Figure 10 shows the comparisons of satellite DCB stabilities between 21st to 30th in January 2011 and 3rd to 12th in May 2015. It can be seen from Figure 10 that the stabilities of satellite DCBs become worse during the solar disturbed period, which corresponds to May 2015, than those during the solar quiet period corresponding to January 2011, especially for PRN 8, 15 and 20 satellites, of which the day scatters are larger than 0.3 ns. Besides, the median value of day scatter increases from 0.08 ns during the solar quiet period to 0.15 ns during the solar disturbed period. There are two possible reasons for the degradation of accuracy, apart from the systematic errors including the mapping error and the ionospheric effective height determination error. One is the overall increase of ionospheric variabilities under stronger and disturbed solar condition, which leads to more errors in local ionospheric modeling and further affects the accuracy of DCB estimates. For example, during the solar disturbed period, the two-degree polynomial may not reproduce the horizontal structure of the local ionospheric VTEC so well as during the solar quiet period. The other reason may be the geographic distribution of the global stations selected. Because there are many more stations in the northern hemisphere, there could be an extra seasonal effect between May and January on the DCB estimates. Figure 11 shows the mean difference and RMS with respect to interpolated CODE GIM over 20 sample stations during the disturbed period. Since there were some stations which did not provide records during the disturbed period, we do not label the missed stations in Figure 11 but leave a blank space for the convenience of comparison with Figure 4. Comparing Figures 4 and 11, it can be seen that the mean differences slightly increase during disturbed period, while RMS values increase significantly during the disturbed period, which means that the time series of local VTEC modeling are more fluctuated during the disturbed period than during the quiet period. This means that the performance of the recursive method is prone to the variabilities of the ionospheric TEC to some degree. Thus, we conclude that the performance of the recursive method can be affected by the solar condition, which is similar to other algorithms such as GIM. More studies are needed to reduce the possible error sources, especially during the disturbed period.

Conclusions
In this paper, we presented a recursive method to estimate satellite DCB and local ionospheric modeling, together with an independent reference station selection procedure. The implementation of the recursive method consists of three parts: first, the local ionospheric VTEC model and DCB of satellites and receivers are estimated based on a recursive filter; second, the receiver DCBs are estimated using the initial satellite DCB estimates from the first step based on the same recursive method; third, the stations which performed worse than the other stations are abandoned from the set of reference stations, according to the difference between two receiver DCB estimates and repeating the procedures above until no station is discarded in the third step. Based on the analyses with respect to the CODE product, the following conclusions are drawn: (1) The accuracy of satellite DCB estimates obtained by the recursive method can reach about 0.10 ns under solar quiet conditions; (2) the influence of stations with bad performances on DCB estimation can be reduced through the independent iterative reference station selection; (3) the accuracy of local ionospheric modeling based on recursive filter can reach less than 2 TECU in the monthly median sense under solar quiet conditions.

Conclusions
In this paper, we presented a recursive method to estimate satellite DCB and local ionospheric modeling, together with an independent reference station selection procedure. The implementation of the recursive method consists of three parts: first, the local ionospheric VTEC model and DCB of satellites and receivers are estimated based on a recursive filter; second, the receiver DCBs are estimated using the initial satellite DCB estimates from the first step based on the same recursive method; third, the stations which performed worse than the other stations are abandoned from the set of reference stations, according to the difference between two receiver DCB estimates and repeating the procedures above until no station is discarded in the third step. Based on the analyses with respect to the CODE product, the following conclusions are drawn: (1) The accuracy of satellite DCB estimates obtained by the recursive method can reach about 0.10 ns under solar quiet conditions; (2) the influence of stations with bad performances on DCB estimation can be reduced through the independent iterative reference station selection; (3) the accuracy of local ionospheric modeling based on recursive filter can reach less than 2 TECU in the monthly median sense under solar quiet conditions.

Conclusions
In this paper, we presented a recursive method to estimate satellite DCB and local ionospheric modeling, together with an independent reference station selection procedure. The implementation of the recursive method consists of three parts: first, the local ionospheric VTEC model and DCB of satellites and receivers are estimated based on a recursive filter; second, the receiver DCBs are estimated using the initial satellite DCB estimates from the first step based on the same recursive method; third, the stations which performed worse than the other stations are abandoned from the set of reference stations, according to the difference between two receiver DCB estimates and repeating the procedures above until no station is discarded in the third step. Based on the analyses with respect to the CODE product, the following conclusions are drawn: (1) The accuracy of satellite DCB estimates obtained by the recursive method can reach about 0.10 ns under solar quiet conditions; (2) the influence of stations with bad performances on DCB estimation can be reduced through the independent iterative reference station selection; (3) the accuracy of local ionospheric modeling based on recursive filter can reach less than 2 TECU in the monthly median sense under solar quiet conditions. It is also necessary to note that the performance of the recursive method is sensitive to solar condition, and more works are needed to further mitigate the dependence of the recursive method on the solar condition. The overall increase of ionospheric variabilities under stronger and more disturbed solar conditions and asymmetry geographic distribution of the global stations selected can be responsible for the degradation of the performance. Considering the low cost and simplicity of implementation, the algorithm has potentials to be implemented in the near real-time DCB estimation and GNSS data quality check in the future [47].
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/6/951/s1, Figure S1: Successive daily changes of geographic distribution of reference stations (shown in green dot) and discarded stations (shown in red dot) during 15th to 30th in January 2011.