#### Abstract

Geographical distribution of global navigation satellite system (GNSS) ground monitoring stations affects the accuracy of satellite orbit, earth rotation parameters (ERP), and real-time 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 ultra-rapid orbit determination and real-time clock offset estimation is conducted. The distribution effects of stations on the products accuracy are analyzed. It shows that the accuracies of GNSS ultra-rapid observed and predicted orbits and real-time 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 time-consuming, 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 real-time clock offset.

#### 1. Introduction

Nowadays, GNSS (global navigation satellite systems) postprocessed and real-time precise satellite orbits and clock offsets, long-term and short-term 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 well-distributed 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 well-distributed 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 5-fold 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 real-time 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) real-time 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 Multi-GNSS 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 multi-GNSS 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 real-time service capabilities are gradually improving. More works are required to optimize the station configuration to quickly produce the ultra-rapid orbit and real-time clock offset high-precision products and to shorten the computation time and to further improve GNSS real-time 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 ultra-rapid orbit and real-time 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 two-dimensional case, the GDOPs reach a minimum at the center of the regular polygon and any combination of them. In the three-dimensional 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 over-size 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 well-distributed stations, improving the accuracy of ultra-rapid orbit determination and real-time clock offset estimation and the calculation of ERP.

#### 2. Station-Geocenter 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 station-geocenter GDOP of the* n* ground stations as an indicator of the uniformity degree of the ground station distribution.

The station-geocenter 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* i*-th station. is the Euclidean distance from the* i*-th 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 station-geocenter 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 three-dimensional 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 station-geocenter 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 low-performance 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 10*n* × 10*n* 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 latitude-longitude grid do not cover equivalent Earth surfaces. In other words, the latitude-longitude grid is not an equal-area grid. From the equator to the poles, the change of latitude-longitude grid area gradually increases. Therefore, the latitude-longitude grid is not evenly divided in the real three-dimensional 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 root-mean-square error in* X*,* Y*, and* Z* directions of the station coordinates, is the probability of the* j*-th point in the* l*-th grid, is the number of points in the grid, is the multipath error value of the* j*-th 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 latitude-longitude 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 high-precision tracking station coordinates in the SINEX (Software INdependent EXchange) format since 1999. The coordinates are obtained through a combination of at least seven AC (analysis-center) 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 root-mean-square 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/data-processing/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 one-step method is used for ultra-rapid orbit determination [25] and the zero-difference method is used for real-time clock offset estimation [26]. The one-step 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 zero-difference 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 real-time clock offset estimation is not made under the condition of accessing a real-time data stream and the support of a related communication chain, but instead uses the ultra-rapid orbit in a file format.

It is noted that the station selection in this paper is only performed in the preprocessing stages of each ultra-rapid and real-time 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 real-time 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 real-time 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 ultra-rapid orbit determination and real-time 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., 24-h 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 ultra-rapid orbit. The ultra-rapid orbit is determined once at GPST 00-h, 06-h, 12-h, and 18-h on each day to generate four 6-h predictions of the ultra-rapid 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 6-h predictions of the ultra-rapid orbit. Using the prediction files and corresponding three-day observation data (day file data; sampling interval of 30 s), the real-time clock offset is solved by employing a square-root information filter algorithm [26, 30]. In this paper, IGS final products are used as a reference for the accuracy statistics of the ultra-rapid orbit and real-time clock offset products of the CUM. It is noted that, for each day, the results of the orbit solution at GPST 00-h and clock offset solution for the whole day are used to derive accuracy statistics.

In this paper, the orbit accuracy is the root-mean-square 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 real-time and final products. The IGS broadcast ephemeris, ultra-rapid (IGU) and the real-time service (RTS) products mainly aim for real-time 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 real-time 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 real-time clock offset estimation, the stochastic model of the clock offset can be described by a first-order Gauss-Markov process. The actual situation shows that the satellite clock offset is mainly represented by high-frequency 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 first-order 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) Single-Day 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, ultra-rapid orbit determination and real-time clock offset estimation are performed, respectively. Experiments conducted for the nine groups are referred to as the first-group, second-group,..., ninth-group 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 third-group, sixth-group, and ninth-group 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 ultra-rapid observed and predicted orbits. Figure 6 shows variations in accuracy (STD) of the real-time 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 real-time 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 ultra-rapid 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 ultra-rapid predicted orbit are, respectively, 59.10 and 31.39 cm, while the differences for the real-time 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 ultra-fast 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 ultra-rapid orbit determination and real-time 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 “ddd-GRID” and “ddd-WSDOP” refer to orbit determination and real-time clock offset estimation with the traditional and probability methods on the “dddth” day, respectively. The “GRID-xx” and “WSDOP-xx” refer to orbit determination and real-time clock offset estimation for “xx” stations with the traditional and probability methods on the 68th day of 2017, respectively.

**(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 real-time 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 ultra-rapid 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 ultra-rapid orbit determination with the LSQ algorithm has a high computational efficiency. Even with 90 stations, ultra-rapid orbit determination with IGS hourly combined observation data can be performed hourly. For the real-time 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 real-time 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 real-time clock offset estimation, the high-frequency 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 real-time clock offset can be updated per 10s. With the corresponding interpolation or forecast methods, the real-time clock can still be updated per second. It can provide a reference for the GPS/BDS fusion real-time clock offset estimation. Because of insufficient stations with BDS real-time data streams at the current stage, it is necessary to combine many stations with IGS real-time data streams to improve the accuracy of BDS real-time 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 ultra-rapid orbit determination and real-time 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 real-time 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 ultra-rapid orbit observation and prediction and the real-time 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, ultra-rapid orbit determination, and real-time 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, real-time 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.