Research Article  Open Access
Random Optimization Algorithm on GNSS Monitoring Stations Selection for UltraRapid Orbit Determination and RealTime Satellite Clock Offset Estimation
Abstract
Geographical distribution of global navigation satellite system (GNSS) ground monitoring stations affects the accuracy of satellite orbit, earth rotation parameters (ERP), and realtime satellite clock offset determination. The geometric dilution of precision (GDOP) is an important metric used to measure the uniformity of the stations distribution. However, it is difficult to find the optimal configuration with the lowest GDOP when taking the 71% ocean limitation into account, because the ground stations are hardly uniformly distributed on the whole of the Earth surface. The station distribution geometry needs to be optimized and besides the stability and observational quality of the stations should also be taken into account. Based on these considerations, a method of configuring global station tracking networks based on grid control probabilities is proposed to generate optimal configurations that approximately have the minimum GDOP. A random optimization algorithm method is proposed to perform the station selection. It is shown that an optimal subset of the total stations can be obtained in limited iterations by assigning selecting probabilities for the global stations and performing a Monte Carlo sampling. By applying the proposed algorithm for observation data of 201 International GNSS Service (IGS) stations for 3 consecutive days, an experiment of ultrarapid orbit determination and realtime clock offset estimation is conducted. The distribution effects of stations on the products accuracy are analyzed. It shows that the accuracies of GNSS ultrarapid observed and predicted orbits and realtime clock offset achieved using the proposed algorithm are higher than those achieved with the traditional method having the drawbacks of lacking evaluation indicators and being timeconsuming, corresponding to the improvements 17.15%, 19.30%, and 31.55%, respectively. Only using 30 stations selected by the proposed method, the accuracies achieved reach 2.01 cm (RMS), 4.93 cm (RMS), and 0.20 ns (STD), respectively. Using 60 stations, the accuracies are 1.47 cm, 3.50 cm, and 0.17 ns, respectively. With the increasing number of stations, the accuracies of the Global Positioning System (GPS) orbit and clock offset improve continuously, but more than 60 stations, the improvement on the orbit determination becomes more gradual, while for more than 30 stations, there is no appreciable increase in the accuracy of the realtime clock offset.
1. Introduction
Nowadays, GNSS (global navigation satellite systems) postprocessed and realtime precise satellite orbits and clock offsets, longterm and shortterm Earth rotation parameters (ERP), interfrequency deviation parameters, and troposphere and ionosphere correction parameters, are produced using the observation data of ground monitoring stations to meet the needs of PNT (positioning navigation and timing) users. Usually we select a subset of welldistributed stations of high quality to obtain the products required by users. Dvorkin and Karutin [1] comprehensively analyzed the effect of the distribution of monitoring stations on the accuracy of GLONASS satellite orbit determination and clock offset estimation. It is shown that 21 welldistributed ground monitoring stations may achieve quadruple coverage of GLONASS satellites (i.e., at least four stations tracking any GLONASS satellite at any time epoch). If the number of stations is increased by 30%, resulting in 5fold coverage, the accuracy of the satellite orbit and clock offset increases only 0.5%. Lou et al. [2] found that the distribution of stations had a nonnegligible effect on realtime clock offset estimation. Li et al. [3] showed that when 25 ground stations were selected to form a triangulated irregular network (TIN), the accuracy of the Global Positioning System (GPS) realtime clock offset was better than 0.3 ns. However, for larger than 25 tracking stations networks, the accuracy was not observed to significantly further increase with the number of stations. Wang et al. [4] point out that the geometric accuracy of ERP estimation is highest when the station coordinates vectors are orthogonal, and the accuracy is proportional to the squared root of the number of stations. Nowadays, there are more than 500 global monitoring stations of the International GNSS Service (IGS). IGS has initiated the MultiGNSS Experiment (MGEX) project (http://mgex.igs.org/) in 2012 to collect and analyze data of emerging new signals and systems [5]. The MGEX network began to provide BDS observations data in 2013. The international GNSS Monitoring and Assessment System (iGMAS) tracking network (http://www.igmas.org/) is developed by China to collect multiGNSS observations and provide satellite orbit and clock, station coordinates and kinds of products for global users [6]. Until now, more than 180 iGMAS/MGEX stations are tracking the BDS2 satellites. As for BDS3 experimental satellites, their B1I and B3I signals are continuously tracked by 26 stations, including 16 iGMAS stations and 10 MGEX stations. The IGS/MGEX/iGMAS realtime service capabilities are gradually improving. More works are required to optimize the station configuration to quickly produce the ultrarapid orbit and realtime clock offset highprecision products and to shorten the computation time and to further improve GNSS realtime product service capabilities.
The grid control method is widely used to obtain uniformly distributed global stations [7]. This method is simply performed by the grid division of the global region, but in fact it is still hard to obtain an optimal solution because of the uniformity of the longitude and latitude grids and the lack of evaluation. Studies [8, 9] pointed out that geometric dilution of precision (GDOP) values being determined by the geometry of the control points and the unknown point reflect the degree of uniformity of the station distribution; i.e., by minimizing the GDOP value, one may obtain a set of stations uniformly distributed. Zhang et al. [7] discussed the effect of the geometric configuration on the accuracy of orbit determination. They experimentally verified that the accuracy of satellite orbit determination is indeed related the orbit dilution of precision (ODOP) of ground stations. Zhang et al. [10] discussed the improvement in the orbit determination due to the accumulation of multiepoch geometric observational information. They concluded that, for an optimal geometric configuration of Beidou navigation satellite, the smaller ODOP indicates the orbit determination of high accuracy. Lou et al. [2] pointed out since the clock parameters were related to the station position, therefore the station distribution would affect the estimation of the satellite clock correction. The experimental results showed that when the station separation was less than 200 km, the correlation between the troposphere parameter and satellite clock offset parameter was strong and the accuracy of satellite clock offset estimation would reach 0.8 ns by using 68 stations. In GNSS navigation and positioning, the GDOP is calculated under the conditions of the control points being the satellites and the unknown point being the receiver. It is well known that the accuracy of receiver position and clock error parameters are related to the GDOP [11, 12]. In the satellite orbit determination, ERP and clock offset estimation, the GDOP of the stations relative to the geocenter is calculated by taking the ground stations as the control points and the geocenter as the unknown point [7, 13, 14]. The above discussion shows optimizing the latter GDOP may improve the accuracy of ultrarapid orbit and realtime clock offset estimation.
As to the importance of the GDOP metric, GDOP minimization has been widely concerned in the past [15]. Recent research suggests that there may be an infinite number of configurations with the lowest GDOP. In the twodimensional case, the GDOPs reach a minimum at the center of the regular polygon and any combination of them. In the threedimensional case, in the center of the regular polyhedron, Cartesian configuration, spiral configuration, Walker configuration, and any combination of them, the GDOPs are minimal [15], and five types of regular polyhedral, including regular tetrahedrons, cubes, regular octahedrons, icosahedrons, and regular dodecahedrons, can be derived. Conical structures have been used for underwater positioning to control the course of ships [16, 17]. Even so, it is still difficult to find all configurations with the lowest GDOP, furthermore, any theoretically optimal configuration cannot be directly applied to the actual situations, such as the distribution constrains of ground stations and the 71% ocean limitation. We also face the problem of selection efficiency, e.g., selecting 60 stations from 500 stations, there are combinations, and the calculation cost is extremely high. This prompts us to propose a probabilistic approach to perform the optimal station selection process.
A stochastic optimization method based on grid control probability proposed in this paper aims for quickly and automatically selecting the list of stations with dominant geometric distribution and station quality. It is believed that this method with potential high efficiency and optimization selection results can provide an optimized solution for the management of the oversize tracking station network (see for example https://magicgnss.gmv.com/StationMonitor/). Currently, the magicGNSS monitoring service of this network is mainly able to monitor each GNSS stations in data level and positioning and receiver environment level. The stations with better global distribution and data quality corresponding the lowest GDOPs or with some problems corresponding the maximum GDOPs can be rapidly obtained using the method. Therefore, a combination of the dynamic monitoring of the GDOP values of the tracking station network and the existing monitoring information in the original service can proved a value reference for the orbit, clock, and ERP solutions, and the geocentric motion.
This paper first applies the grid method to the candidate stations to allocate the initial probabilities, and then, considering the quality factors of the stations, probabilities are reassigned to the candidate stations. A comprehensive indicator for measuring the uniformity of distribution of ground stations and the observation quality is introduced to justify the selected subset of stations. It is shown that the proposed approach can quickly and automatically select high quality and welldistributed stations, improving the accuracy of ultrarapid orbit determination and realtime clock offset estimation and the calculation of ERP.
2. StationGeocenter DOP Measuring the Uniformity of the Stations
When the stations as known points are evenly distributed on the Earth surface, the GDOP at the geocenter can get the minimum [8]. This indicates that by the stationgeocenter GDOP of the n ground stations as an indicator of the uniformity degree of the ground station distribution.
The stationgeocenter positioning configuration can be recorded as . Here, denotes stations, is geocenter, and m is the dimension of the configuration (m = 3) [10]. For the configuration comprising n ground stations, if and only if [15, 18, 19]wherewhere is the line of sight (LOS) from the geocenter to the ith station. is the Euclidean distance from the ith station to the geocenter x, j is the coordinate component. I is the identity matrix, the GDOP can reach the theoretical minimum [15]Because the GDOP is usually a concept of positioning and cannot be directly applicable to the subject under our discussion. We can define the stationgeocenter DOP asto measure the uniformity of the stations, where tr denotes the trace of the matrix.
Considering the GDOP minimization condition expressed in (1) and (3), the SDOP minimization indicates that [19]andwhere , n; k, s =1, 2 and 3; a_{i,k} or a_{i,k} is the element of the vector a_{i. }The condition (5) results in balanced station geometry to improve the reliability of the GNSS products to compensate the systematic errors. The orthogonality and uniformity conditions in (6) indicate that the total observational information will be optimally allocated to the eight quadrants in threedimensional space; thus, the optimality of the ground stations can be achieved at least for the ERP determination, c.f. Wang et al. [4]
As is well known, the geometry of the station distribution is only one of impact factors affecting the GNSS products, e.g., some of stations improve the geometric configuration obviously, but their quality may be poor, leading to bad GNSS products. For this, we further introduce the following weighted stationgeocenter DOP aswhere is the weight matrix about the stations to reflect the station stability and observational quality, where s_{1} and s_{2} are, respectively, the medians of all station variances and multipath errors. b and d are the empirical values, respectively, 0.7 and 0.3 in this paper. This method not only considers the quality and geometry of the stations but also achieves a good combination of the two above factors, ensuring the orbit determination and clock offset estimation results to be optimal. The weight P considers the theory of robust estimation [20]. When is greater than 1, the weight is directly set to 1. While it is less than 1, it is used according to its actual value. Some points improve the geometric configuration obviously, but there quality is poor, leading the weight of the stations being given too small. Some points improve the geometric configuration badly, but their quality is good, leading the weight of the stations being given too large. The robust process can avoid the above conditions leading the overall WGDOP value being too large, losing the value of the stations in the orbit determination, and thereby reducing the effectiveness of the indicator.
In general, the WSDOP would lead to an overall indicator by combining the station geometry with the station quality factors, thereby an effectiveness indicator.
3. Random Selection Algorithm Based on WSDOP
3.1. Random Optimization Model
The station selection is performed to solve two problems: (1) model incompleteness, such as that due to system error (i.e., station stability or observation quality) and (2) compromise between the calculation cost and the accuracy (i.e., the fast positioning of lowperformance terminals or fast automatic selection when a certain accuracy is guaranteed). However, the above SDOP or WSDOP minimization may face with various limitations and the NP problem, i.e., optimally selecting P points from N points. For this, we propose a probabilistic approach to search the optimal station configuration by minimizing WSDOP; that is,where n denotes the number of all candidate stations. p denotes the number of stations selected. c is station. is a probability space defined aswhere . Because the station is random, the SDOP or WSDOP becomes a random variable. Next, we will show how to perform the random experiment to obtain the lowest WSDOP by allocating the selection probabilities to each candidate stations in the selection algorithm design. It can be noted that because a constant can be thought of as a random variable with zero variance, the traditional method is a special case of the random experiment in this paper.
3.2. Random Selection Algorithm Design
(1) Probability Allocation Approach. Presently, the grid method is usually used to select stations. The basic idea is to divide the global region using a suitable latitude and longitude grid according to the number of selected stations and divides multiple discrete stations into different areas. When the region is divided into a 10n × 10n grid (), for example, the grid number num can be calculated asIn the above formula, denotes rounding up and n = 2, 3, 4 can be used for grid station selection for up to 150 stations [21]. It should be noted that a unit of longitude at a certain latitude is not equivalent in distance to the same unit of longitude at a different latitude; hence cells in a homogeneous latitudelongitude grid do not cover equivalent Earth surfaces. In other words, the latitudelongitude grid is not an equalarea grid. From the equator to the poles, the change of latitudelongitude grid area gradually increases. Therefore, the latitudelongitude grid is not evenly divided in the real threedimensional space, resulting in the fact that geometric strength of the split graphs is theoretically not optimal [22].
Owing to the extremely uneven global distribution of the stations, the number of stations is different in each cell. At this time, it is generally necessary to comprehensively consider the position and quality of each station and to select the ideal orbital station configuration in each cell. To improve the stochastic optimization convergence ability and to avoid the selected stations being concentrated in a certain area, it is necessary to control the allocation probabilities of stations by the grid. On the basis of (11), the global stations are evenly partitioned with latitude and longitude spacing bywhere the unit of JW is degrees and n is the number of stations, and the quantity and quality of the stations in a cell are counted.
The method mainly considers the distribution of stations, and it is difficult to weigh the information of station stability and station observation quality at the same time. There are also many artificial control factors in the selection process. The factors usually refer to artificially determining the appropriate number of grid division by multiple grid scaling, and artificially selecting one station from multiple stations in each grid because of lacking corresponding evaluation indicators. Addressing these drawbacks, we propose a stochastic combinatorial optimization algorithm to obtain the optimal solution from the huge number of candidate configurations for all combinations. On the basis of the above problems, a method of station screening under the grid control probability is used, and the idea of probability statistics is introduced into the station selection. That is to say, on the basis of the grid method and by considering the data quality, stability, and geographical distribution of the station, each station is assigned a certain probability. A station is assigned a high probability when the stability and observation quality of the station are high. The probability allocation formula readswhere is the rootmeansquare error in X, Y, and Z directions of the station coordinates, is the probability of the jth point in the lth grid, is the number of points in the grid, is the multipath error value of the jth point, and b and d are empirical values, respectively, 0.7 and 0.3 in this paper.
(2) Monte Carlo Experiments and Best Solution Searching. According to the principle of Monte Carlo random sampling, a corresponding probability is first assigned to each station:where num_block is the number of cells containing a or more stations. That is, after the global region being divided by latitudelongitude grids according to formula (12), the grids containing no stations are removed, and the number of the remaining grids is the num_block. F_{j} is the probability of station j.
Taking all stations as the experimental population in each random sampling, n stations are randomly selected from all the candidate stations S (i.e., a list of stations is randomly selected): where F is a matrix containing probabilities of all stations, randsrc denotes an algorithm for selecting n stations from the total stations S according to the probability F, and select_list denotes a station list. At the same time, the WSDOP is calculated in each sampling. In order to obtain a better station list, it is needed to perform multisampling, such as 100,000 times. Finally, by evaluating the size of WSDOP of each sample, we can obtain a solution with the lowest WSDOP among the huge samples.
4. Experimental Analysis
4.1. Algorithm Implementation and Process Design
Experimental data for a total of 201 ground stations are downloaded from the IGS website. IGS has provided a highprecision tracking station coordinates in the SINEX (Software INdependent EXchange) format since 1999. The coordinates are obtained through a combination of at least seven AC (analysiscenter) products [23]. For the stations without precision coordinates, such as the iGMAS stations and the regional continuously operating reference stations (CORS), their precision coordinates can be obtained using precise point positioning (PPP) technique. Therefore, the precise coordinates and their rootmeansquare error for the stations in this paper are obtained from the IGS SINEX file. The multipath error values are calculated with Translation Editing and Quality Checking (TEQC) which is an internationally recognized and widely used GNSS observation data preprocessing software (https://www.unavco.org/software/dataprocessing/teqc/teqc.html). The data quality checking (QC) is the core function of TEQC. But it can only analyze the GPS/GLONASS observation data with RINEX2.XX or earlier formats. It cannot analyze the observation data of BDS/GALILEO. For the shortage of TEQC, the quality analysis software being developed dy BDS Data Processing Center at the China University of Mining and Technology (CUM) can process multisystem observation data with RINEX2.XX and RINEX3.XX formats, including BDS3 experiment satellite data [24]. After randomly selecting stations according to the grid control probability (with the number of samples set to 100,000), a onestep method is used for ultrarapid orbit determination [25] and the zerodifference method is used for realtime clock offset estimation [26]. The onestep method based on nondifferential measurement can effectively reduce the influence of system error, especially in multisystem fusion orbit determination [25, 27]. Ambiguity parameters can be preserved using the zerodifference method, which provides the possibility of subsequent double difference ambiguity fixation and satellite hard delay estimation [28, 29]. Figure 1 shows the flow of the station selection algorithm and the orbit and clock offset solution process.
Figure 1 shows that (1) the list of stations corresponding to the smallest WSDOP value in the station selection process is the final station selection result and (2) the realtime clock offset estimation is not made under the condition of accessing a realtime data stream and the support of a related communication chain, but instead uses the ultrarapid orbit in a file format.
It is noted that the station selection in this paper is only performed in the preprocessing stages of each ultrarapid and realtime clock offset solutions rather than the processing stages. The weight P is designed with the idea of the robust estimation to ensure a more stable selection results. With the continuous development of artificial intelligence (AI) and IGS realtime services, and the increasing of the ground stations, such as the global iGMAS monitoring stations and the CORS stations, the stations selection method evolution is needed to be further considered with a more operational strategy. For example, the neural network or the ant colony algorithm can be used to further improve the search efficiency of the global optimal solution of the method. The number of sampling of Monte Carlo random experiments can be also reduced under the condition of ensuring a certain precision, such as 2000 times. The validity of frequent station selection (i.e., the dynamic selection method) in the realtime or the final precision product processing, such as selecting per 10 minutes, should be further studied.
4.2. Data Processing and Experimental Analysis
The experiment conducted in this paper mainly investigates the effects of the station distribution and number on the accuracy of ultrarapid orbit determination and realtime clock offset estimation using the WSDOP value. After optimizing the station selection, using the orbit and clock software of CUM, IGS observation data (i.e., 24h observation arc data merged with hourly observation files) for 3 consecutive days, from the 68th day to the 70th day of 2017, are used for the determination and prediction of the ultrarapid orbit. The ultrarapid orbit is determined once at GPST 00h, 06h, 12h, and 18h on each day to generate four 6h predictions of the ultrarapid orbit. When 90 stations are used in orbit determination on each day, the orbital prediction file of the day is generated by concatenating the four 6h predictions of the ultrarapid orbit. Using the prediction files and corresponding threeday observation data (day file data; sampling interval of 30 s), the realtime clock offset is solved by employing a squareroot information filter algorithm [26, 30]. In this paper, IGS final products are used as a reference for the accuracy statistics of the ultrarapid orbit and realtime clock offset products of the CUM. It is noted that, for each day, the results of the orbit solution at GPST 00h and clock offset solution for the whole day are used to derive accuracy statistics.
In this paper, the orbit accuracy is the rootmeansquare error (RMSE), while the clock accuracy is the mean Standard Deviations of Errors (STDE). The twice difference comparison method calculated in two steps is used to obtain the STDE [31, 32]. The flow of the method is that: first, the respective single difference values of the two products are obtained. That is, a satellite is selected as the base satellite, and the single difference values are obtained from base satellite and other satellites, thereby eliminating the difference generated by the effect of the difference in the reference clock on clock offset. Second, the second difference is obtained by the two products of the single difference, which can effectively reflect the coincidence degree of clock offset estimation with different software and strategies. IGS can provide realtime and final products. The IGS broadcast ephemeris, ultrarapid (IGU) and the realtime service (RTS) products mainly aim for realtime applications, while the IGS rapid (IGR) and final products for postmission applications. The nominal accuracies of broadcast orbit and clock are reported as 1.00 m and 2.50 ns, respectively. The nominal accuracies of IGU predicted orbit and clock are 5.00 cm and 1.50 ns, respectively. The nominal accuracies of IGU observed orbits and clocks are 3.00 cm and 3.00 ns, respectively (http://www.igs.org/products). Zhang et al. [33] performed an evaluation and analysis of RTS orbits and clocks products from different IGS analysis centers. The results showed that the RMSE of satellite orbit ranged from 3.80 cm to 7.50 cm for different RTS products. While the mean STDE of satellite clocks ranged from 0.06 ns to 0.19 ns.
The calculation strategy of the orbit determination and realtime clock offset estimation is presented in Table 1. In the strategy, the stochastic function model of the process noise of the unknown parameters is important. It reflects the accuracy of the motion model to a certain extent, and directly affects the contribution of kinematic information to the parameters solution. A suitable stochastic model is established on the analysis of the statistical characteristics of the process noise [30]. Taking an example of the clock offset estimation and modeling, the stochastic function model is analyzed as follows: (1) In the process of the realtime clock offset estimation, the stochastic model of the clock offset can be described by a firstorder GaussMarkov process. The actual situation shows that the satellite clock offset is mainly represented by highfrequency signals, and its process noise can be described by white noise model [25, 34, 35]. The size of the process noise can be set by the analysis of the frequency stability of the satellite clock offset with the phase noise spectral density in the frequency domain or the Allan variance in the time domain. (2) In the process of satellite clock offset modeling, Wang et al. [36] analyzed the clock offset system, periodic and random characteristics. The clock was modeled by firstorder stochastic differential equations and its optimal estimation was realized by Kalman filtering algorithm. The state parameters of the basic model of the clock model include the clock offset, relative frequency deviation, and frequency drift, and the corresponding random noises are frequency white noise, frequency random walk noise and frequency random running noise. Considering the periodic components, phase and frequency flicker noises, the basic model was further extended.

(1) SingleDay Station Selection Experiment of the Grid Control Probability. Nine groups of stations are selected from the stations on the 68th day, with the number of stations in each group being . The WSDOP values in ascending order and the corresponding SDOP values according to (4) for each group of stations are obtained. After selecting the lists of stations corresponding to the first 10 and last 10 values of each group of WSDOP values, ultrarapid orbit determination and realtime clock offset estimation are performed, respectively. Experiments conducted for the nine groups are referred to as the firstgroup, secondgroup,..., ninthgroup experiments. Figure 2 shows the variations of the actual minimum WSDOP value, minimum SDOP, and theoretical minimum GDOP (Equation (3)) after optimized selection of the nine groups stations on the 68th day, with an increase in the number of stations. Figure 3 shows the variations in WSDOP and SDOP values for each group’s experiment. Figure 4 shows the frequency histograms and probability density functions of WSDOP and SDOP values for the thirdgroup, sixthgroup, and ninthgroup experiments.
Figure 2 shows that, in general, the WSDOP and theoretical GDOP values gradually approach with an increase in the number of stations. After the number of stations reaches a certain number, the curves of the two values eventually coincide. As the number of stations increases, the WSDOP, SDOP, and theoretical GDOP gradually decrease; i.e., the geometry of the stations gradually improves. When there are few stations, especially when there are fewer than 60 stations, the change rate in the GDOP value is large, and the amplitude of the variation in the WSDOP value is larger than the amplitudes of the variations in the theoretical GDOP and SDOP values. This shows that the geometric configuration significantly improves as the number of stations increases within a certain range. When there are more than 60 stations, the change in the three kinds of DOP values is almost negligible, and their rate of change tends to zero. In this case, therefore, the improvement in the geometry of the stations gradually slows as the number of stations increases.
Figure 3 shows that the SDOP value corresponding to the WSDOP minimum is not necessarily the smallest, but the overall varying trend of the two values is consistent. This indicates that the WSDOP value simultaneously considers the geometric distribution and quality of the stations. Additionally, the fewer the stations, the greater the change in WSDOP and SDOP values in each group’s experiment, indicating the more obvious advantage of using the indicator WSDOP in station selection.
Figure 4 shows that, for the experiments of the three groups, the distributions of WSDOP and SDOP values obtained in this paper can be modeled with normal distributions, indicating the random experiments are successfully performed.
Twenty experiments of orbit and clock offset determination are conducted for each group above. Figure 5 shows variations in accuracy (RMS) of the ultrarapid observed and predicted orbits. Figure 6 shows variations in accuracy (STD) of the realtime clock offset.
Figures 3, 5, and 6 provide three findings. (1) When the number of stations is the same, the accuracy of orbit determination and realtime clock offset estimation is better in the first 10 experiments than in the last 10 experiments, and the WSDOP value indicates the accuracy of the orbit and clock offset to some extent; i.e., for the same number of stations, the accuracy of the orbit and clock offset improves as the WSDOP value decreases. This accuracy improvement is more obvious when there are few stations (i.e., 30 stations or fewer). When there are 10 stations, the difference between the highest accuracy in the first 10 experiments and the lowest accuracy in the last 10 experiments for the ultrarapid observed orbit is 23.60 cm while the difference between the average accuracy in the first 10 experiments and the average accuracy in the last 10 experiments is 9.24 cm; the differences for the ultrarapid predicted orbit are, respectively, 59.10 and 31.39 cm, while the differences for the realtime clock offset are, respectively, 0.81 and 0.21 ns. When there are 30 stations, the accuracy differences for the observed orbit are, respectively, 6.38 and 1.99 cm, the differences for the orbit predicted orbit are 11.10 and 6.37 cm, and the differences for the clock offset are 0.40 and 0.20 ns. When there are 90 stations, the differences for the observed and predicted orbits and the clock offset are, respectively, 0.67 and 0.22 cm, 1.40 and 0.96 cm, and 0.16 and 0.00 ns. The results show that when there are more than a certain number of stations, the effect of the geometric configuration on the accuracy of the orbit and clock offset gradually weakens. (2) The fluctuation in the calculation results of the first 10 experiments is small while the fluctuation in the last 10 runs is larger. It is concluded that, for the same number of stations, the result is more stable when the WSDOP value is smaller.
Figures 2, 3, 4, 5, and 6 provide four findings. (1) An optimal subset of the total stations can be obtained in limited iterations by assigning selecting probabilities for the global stations and performing the Monte Carlo sampling. (2) The WSDOP value well reflects the accuracy of the orbit and clock offset and the accuracy has a negative correlation with the WSDOP value. (3) The accuracies of the ultrafast observed and predicted orbits are consistent; i.e., the predicted orbit is good when the accuracy of the observed orbit is good. (4) Improvements in the accuracy of the orbit and clock offset are consistent; i.e., the clock offset improves where the accuracy of the orbit improves.
(2) Multiday Station Selection Experiments of Grid Control Probability and Traditional Grid Methods. Using the probability method and traditional method (Equation (12)), nine groups of experimental station selections are executed, with the number of stations in each group being 10, 20,…, 90. The ultrarapid orbit determination and realtime clock offset estimation are then performed with data for 3 consecutive days. Figures 7(a), 7(b), and 7(c), respectively, presents the station selection results obtained using the two methods for 30, 60, and 90 stations on the 68th day. The names and longitudes and latitudes of the station selection lists are shown in the Appendix. Figure 8 shows the accuracy statistics of the nine groups of experiments for 3 consecutive days using the two methods. Table 2 presents the corresponding average accuracy and calculation time statistics of the 3 consecutive days. In the following, the “dddGRID” and “dddWSDOP” refer to orbit determination and realtime clock offset estimation with the traditional and probability methods on the “dddth” day, respectively. The “GRIDxx” and “WSDOPxx” refer to orbit determination and realtime clock offset estimation for “xx” stations with the traditional and probability methods on the 68th day of 2017, respectively.
 
Note: the station selection is performed with a conventional computer while the orbit and clock calculations are made with a graphics workstation. 
(a)
(b)
(c)
(a)
(b)
(c)
(1) Figure 8 and Table 2 show that the accuracy of the stochastically optimized station selection method is higher than that of the traditional grid method. The accuracies of the orbit and clock offset have the highest rates of increase when there are 10 stations. When there are 50 stations, the average accuracies of the orbit and clock offset achieved with the grid method are, respectively, 1.87 cm, 4.67 cm, and 0.26 ns while the accuracies achieved with the probability method are better by 7.14%, 11.43%, and 27.85%, respectively. When there are 90 stations, the accuracies achieved with the grid method are, respectively, 1.27 cm, 3.53 cm, and 0.23 ns while the accuracies achieved with the probability method are better by 9.75%, 11.32%, and 23.19%, respectively. At the same time, the advantages of the probability method are more obvious for fewer stations (e.g., 20 stations or fewer). The positive effect of the accuracy improvement is better for the clock offset than for the orbit. As a whole, the accuracies achieved with the probability method can reach the same level as IGS RTS products [33].
(2) Table 2 shows that, on the whole, the accuracies of the GPS orbit and clock offset increase gradually as the number of stations increases. When the probability method is applied to the selection of 30 stations, the accuracies of the observed and predicted orbits and clock offset are, respectively, 2.01 cm, 4.93 cm, and 0.20 ns. When there are fewer than 60 stations, the accuracies of the orbit and clock offset appreciably improve with an increase in the number of stations. When there are more than 60 stations, the accuracy of the orbit does not improve greatly. When there are more than 30 stations, the improvement in the accuracy of the realtime clock offset is no longer obvious. When there are more than 60 stations, there is a slight decrease in the accuracy of the orbit and the clock offset with an increase in the number of stations.
(3) Figure 7 and Table 2 show that the grid control probability method considers the advantages of the uniform distribution of the traditional grid method and also the quality of the station observation data and station stability, thus improving accuracy.
(4) Table 2 shows that, with the probability method, the maximum calculation time of station selection does not exceed 14.14 min. The calculation time for the ultrarapid orbit and clock offset increases gradually with the number of stations, and, in particular, the calculation time of the clock offset sharply increases. The calculation time of the orbit with 90 stations is 1.48 times that with 60 stations. The calculation time of the clock offset with 90 stations is 29.04 times that with 30 stations. It is therefore necessary to select an appropriate number of stations when meeting the accuracy requirements of the orbit and clock offset. It can be seen from the computation time of the orbit determination and clock offset estimation that the ultrarapid orbit determination with the LSQ algorithm has a high computational efficiency. Even with 90 stations, ultrarapid orbit determination with IGS hourly combined observation data can be performed hourly. For the realtime clock offset estimation with the SRIF algorithm, the computation time of each epoch solution with 50 stations is 1.17 s which is obtained by the total time 53.43min being divided by the total epoch 2880. It indicates that the realtime clock offset estimation in this case cannot achieve the clock offset update per second. Considering the network delay and other factors affecting the timeliness of the realtime clock offset estimation, the highfrequency clock offset estimation with no more than 40 stations are recommended. At the same time, when the number of stations is 90, the computation time of each epoch clock offset solution is 7.50s. It indicates that the realtime clock offset can be updated per 10s. With the corresponding interpolation or forecast methods, the realtime clock can still be updated per second. It can provide a reference for the GPS/BDS fusion realtime clock offset estimation. Because of insufficient stations with BDS realtime data streams at the current stage, it is necessary to combine many stations with IGS realtime data streams to improve the accuracy of BDS realtime clock offset [26].
5. Conclusion
A random optimization algorithm for the selection of GNSS global monitoring stations based on the selection judgment index WSDOP was proposed. After analyzing optimized ground station distribution obtained by the proposed probabilistic algorithm, it is shown that the proposed algorithm quickly and automatically produces a desirable subset of total stations with regard to the geometric distribution factor and the quality factors about the station observations. By applying this algorithm to the calculation of ultrarapid orbit determination and realtime clock offset, we find that the accuracies of the orbit and clock offset for 3 consecutive days achieved are better than that of the traditional grid method. The conclusions are drawn as follows.
(1) In general, irrespective of whether the number of stations is the same or different, while the geometric configuration of the stations is improved, the accuracy of the orbit and realtime clock offset are both improved with the decreasing WSDOP values.
(2) There is an obvious improvement in the station geometric configuration with an increase in the number of stations. The improvement in the orbital geometric structure information is gradual when there are more than 60 stations. When using 30 stations selected by the proposed method, the accuracies of the GPS ultrarapid orbit observation and prediction and the realtime clock offset can reach 2.01 cm, 4.93 cm, and 0.20 ns, respectively. When using 30 stations, the accuracies are, respectively, 1.47 cm, 3.5 cm, and 0.17 ns. For more than 60 stations, the accuracy improvement of the orbit determination will be limited; meanwhile, 30 stations may achieve an extreme accuracy of the clock estimation.
(3) The GNSS product accuracies from the proposed algorithm are overall higher than those of the traditional grid method, with the observation and prediction orbit and clock offset accuracies improving by 4.12–59.16%, 7.08–56.33%, and 17.65–58.43% over the latter method with the average increasing rates of 17.15%, 19.30%, and 31.55%. When few stations are selected (e.g., 20 stations or fewer), the improvements are more remarkable both in the orbit and clock offset.
(4) The calculation times of station selection, ultrarapid orbit determination, and realtime clock offset estimation increase gradually with the number of stations; in particular, the calculation time of the clock offset sharply increases. When the number of station selection is 90 stations are selected, the respective calculation times are 14.15, 16.82, and 358.04 min.
Combined with the neural network method or ant colony algorithm, the probability method presented in this paper is expected to be further improved. Considering that the monitoring stations of the Beidou system are unevenly and less distributed presently, the method can provide an effective approach for orbit determination, realtime clock offset estimation, ERP solutions, and construction of global stations of the Beidou system.
Appendix
See Table 3.

Data Availability
The data used to support the findings of this study are available from the first or the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
Qianxin Wang, Shuqiang Xue, and Xu Yang conceived and defined the research scheme. Xu Yang and Shuqiang Xue verified the feasibility of the method and implemented the software algorithm; Xu Yang and Qianxin Wang checked the data processing results and wrote the manuscript; Qianxin Wang and Shuqiang Xue helped to revise the manuscript and modified some figures and tables.
Acknowledgments
This work was supported by “the National Science and Technology Basic Work of China (no. 2015FY310200)”, “the State Key Program of National Natural Science Foundation of China (no. 41730109)”, “the National Natural Science Foundation of China (nos. 41404033 and 41874039)”, and “the Jiangsu Dual Creative Teams Program Project Awarded in 2017”; and IGS and iGMAS are acknowledged for providing data.
References
 V. V. Dvorkin and S. N. Karutin, “Optimization of the global network of tracking stations to provide GLONASS users with precision navigation and timing service,” Gyroscopy and Navigation, vol. 4, no. 4, pp. 181–187, 2013. View at: Publisher Site  Google Scholar
 Y. Lou, X. Dai, and W. Song, “Research on the influence of stations' distance in highaccuracy GPS satellite clock offset estimation,” Geomatics & Information Science of Wuhan University, vol. 36, no. 4, pp. 397–406, 2011. View at: Google Scholar
 H. Li, Y. Li, T. Dai, L. Wang, and J. Wang, “Realtime GPS satellite clock offset determination based on the TIN global stationselecting method,” Journal of Geomatics Science & Technology, vol. 34, no. 5, pp. 441–444, 2017. View at: Google Scholar
 Q. Wang, Y. Dang, and T. Xu, “The method of earth rotation parameter determination using GNSS observations and precision analysis,” Lecture Notes in Electrical Engineering, vol. 243, pp. 247–256, 2013. View at: Publisher Site  Google Scholar
 O. Montenbruck, P. Steigenberger, L. Prange et al., “The multiGNSS experiment (MGEX) of the international GNSS service (IGS) – achievements, prospects and challenges,” Advances in Space Research, vol. 59, no. 7, pp. 1671–1697, 2017. View at: Publisher Site  Google Scholar
 H. Cai, G. Chen, W. Jiao, K. Chen, T. Xu, and H. Wang, “An initial analysis and assessment on final products of iGMAS,” in China Satellite Navigation Conference (CSNC) 2016 Proceedings, vol. 3 of Lecture Notes in Electrical Engineering, pp. 515–527, Springer, Singapore, 2016. View at: Publisher Site  Google Scholar
 R. Zhang, Y. Yang, G. Huang, L. Wang, W. Qu, and H. Zhao, “Optimal solution of dynamic parameters of BDS satellite and optimal distribution of orbit determination tracking station,” Journal of Geodesy and Geodynamics, vol. 36, no. 3, pp. 216–220, 2016. View at: Google Scholar
 S. Xue and Y. Yang, “Nested cones for singlepointpositioning configuration with minimal GDOP,” Geomatics and Information Science of Wuhan University, vol. 39, no. 11, pp. 1369–1374, 2014. View at: Google Scholar
 J. Li, Z. Li, W. Zhou, and S. Si, “Study on the minimum of GDOP in satellite navigation and its applications,” Acta Geodaetica et Cartographica Sinica, vol. 40, pp. 85–88, 2011. View at: Google Scholar
 L. Zhang, Y. Dang, Y. Cheng, S. Xue, S. Gu, and D. Han, “Analysis and optimization of BDS GEO/IGSO/MEO ground monitoring stations configuration for determining GNSS orbit,” Acta Geodaetica et Cartographica Sinica, vol. 45, no. S2, pp. 82–92, 2016. View at: Google Scholar
 Y. Yang, “Progress, contribution and challenges of compass/beidou satellite navigation system,” Acta Geodaetica et Cartographica Sinica, vol. 39, no. 1, pp. 1–6, 2010. View at: Google Scholar
 Y. X. Yang, J. L. Li, J. Y. Xu, J. Tang, H. R. Guo, and H. B. He, “Contribution of the compass satellite navigation system to global PNT users,” Chinese Science Bulletin, vol. 56, no. 26, pp. 2813–2819, 2011. View at: Publisher Site  Google Scholar
 A. M. ElNaggar, “New method of GPS orbit determination from GCPS network for the purpose of DOP calculations,” Alexandria Engineering Journal, vol. 51, no. 2, pp. 129–136, 2012. View at: Publisher Site  Google Scholar
 X. Peng, Research on Optimization Method of Orbit Determination Station Based on DOP Value, China University of Mining and Technology, 2016.
 S. Xue and Y. Yang, “Positioning configurations with the lowest GDOP and their classification,” Journal of Geodesy, vol. 89, no. 1, pp. 49–71, 2015. View at: Publisher Site  Google Scholar
 V. Ballu, M.N. Bouin, S. Calmant et al., “Absolute seafloor vertical positioning using combined pressure gauge and kinematic GPS data,” Journal of Geodesy, vol. 84, no. 1, pp. 65–77, 2010. View at: Publisher Site  Google Scholar
 J. Zhao, Y. Zou, H. Zhang, Y. Wu, and S. Fang, “A new method for absolute datum transfer in seafloor control network measurement,” Journal of Marine Science and Technology, vol. 21, no. 2, pp. 216–226, 2016. View at: Publisher Site  Google Scholar
 S. Xue, Y. Yang, and Y. Dang, “A closedform of newton method for solving overdetermined pseudodistance equations,” Journal of Geodesy, vol. 88, no. 5, pp. 441–448, 2014. View at: Publisher Site  Google Scholar
 S. Xue and Y. Yang, “Understanding GDOP minimization in GNSS positioning: infinite solutions, finite solutions and no solution,” Advances in Space Research, vol. 59, no. 3, pp. 775–785, 2017. View at: Publisher Site  Google Scholar
 Y. Yang, “Analysis and analysis of geocentric motion by the use of comprehensive power,” China Surveying and Mapping Society Science and Technology Information Network Branch, vol. 7, 2002. View at: Google Scholar
 T. Dai, J. Li, J. Zhao, P. Pang, and Y. Wei, “Study of global stationselection and application in precise orbit determination based on TIN NET,” Journal of Geodesy Geodynamics, vol. 37, no. 1, pp. 77–80, 2017. View at: Google Scholar
 S. Xue, Research on Geodetic Observation Optimization Theory and Methods, Chang'an University, 2018.
 M. Li and T.H. Xu, “Research on the combination of IGS analysiscenter solution for station coordinates and ERPs,” Lecture Notes in Electrical Engineering, vol. 305, no. 3, pp. 15–29, 2014. View at: Publisher Site  Google Scholar
 Y. He, Q. Wang, Z. Wang, and Y. Mao, “Quality analysis of observation data of beidou3 experimental satellites,” in China Satellite Navigation Conference (CSNC) 2018 Proceedings, vol. 498 of Lecture Notes in Electrical Engineering, pp. 275–294, Springer, Singapore, 2018. View at: Publisher Site  Google Scholar
 M. Li, C. Shi, Q. Zhao, and J. Liu, “MultiGNSS precision orbit determination,” Acta Geodaetica et Cartographica Sinica, vol. 40, no. 5, pp. 26–30, 2011. View at: Google Scholar
 Q. Zhao, Z. Dai, G. Wang, X. Li, and J. Liu, “Realtime precise BDS clock estimation with the undifferenced observation,” Geomatics and Information Science of Wuhan University, vol. 41, no. 5, pp. 686–691, 2016. View at: Google Scholar
 P. Pang, J. Li, T. Dai, and Y. Li, “Data quality and accuracy analysis of BDS/GPS nondifferential onestep multiorbit determination,” Journal of Geomatics Science and Technology, vol. 34, no. 6, pp. 593–601, 2017. View at: Google Scholar
 A. Hauschild and O. Montenbruck, “Kalmanfilterbased GPS clock estimation for near realtime positioning,” GPS Solutions, vol. 13, no. 3, pp. 173–182, 2009. View at: Publisher Site  Google Scholar
 X. Li, Y. Xu, and L. Wang, “Undifferenced precise satellite clock error estimation and precision analysis,” Geomatics and Information Science of Wuhan University, vol. 35, no. 6, pp. 661–664, 2010. View at: Google Scholar
 X. Dai, RealTime Precise GNSS Satellite Orbit Determination Using the SRIF Method: Theory and Implemencation, Wuhan University, 2016.
 Y. Lou, C. Shi, X. Zhou, and S. Ye, “Realization and analysis of GPS precise clock products,” Geomatics & Information Science of Wuhan University, vol. 34, no. 1, pp. 88–91, 2009. View at: Google Scholar
 L. Chen, W. Song, W. Yi, C. Shi, Y. Lou, and H. Guo, “Research on a method of realtime combination of precise GPS clock corrections,” GPS Solutions, vol. 21, no. 1, pp. 187–195, 2017. View at: Publisher Site  Google Scholar
 L. Zhang, H. Yang, Y. Gao, Y. Yao, and C. Xu, “Evaluation and analysis of realtime precise orbits and clocks products from different IGS analysis centers,” Advances in Space Research, vol. 61, no. 12, pp. 2942–2954, 2018. View at: Publisher Site  Google Scholar
 M. Li, Researeh on MultiGNSS Precise Orbit Determination Theory and Application, Wuhan University, 2011.
 G. Dai, Research on the GNSS Satellite Clock Offset RealTime Determination and The RealTime Precise Point Positioning, PLA Information Engineering University, 2017.
 B. Wang, Analysis of BDS Satellite Clock in Orbit, Modelling and Its Prediction Research, Wuhan University, 2016.
Copyright
Copyright © 2019 Xu Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.