A Method for Rapidly Determining the Optimal Distribution Locations of GNSS Stations for Orbit and ERP Measurement Based on Map Grid Zooming and Genetic Algorithm

: Designing the optimal distribution of Global Navigation Satellite System (GNSS) ground stations is crucial for determining the satellite orbit, satellite clock and Earth Rotation Parameters (ERP) at a desired precision using a limited number of stations. In this work, a new criterion for the optimal GNSS station distribution for orbit and ERP determination is proposed, named the minimum Orbit and ERP Dilution of Precision Factor (OEDOP) criterion. To quickly identify the specific station locations for the optimal station distribution on a map, a method for the rapid determination of the selected station locations is developed, which is based on the map grid zooming and heuristic technique. Using the minimum OEDOP criterion and the proposed method for the rapid determination of optimal station locations, an optimal or near-optimal station distribution scheme for 17 newly built BeiDou Navigation Satellite System (BDS) global tracking stations is suggested. To verify the proposed criterion and method, real GNSS data are processed. The results show that the minimum OEDOP criterion is valid, as the smaller the value of OEDOP, the better the precision of the satellite orbit and ERP determination. Relative to the exhaustive method, the proposed method significantly improves the computational efficiency of the optimal station location determination. In the case of 3 newly built stations, the computational efficiency of the proposed method is 35 times greater than that of the exhaustive method. As the number of stations increases, the improvement in the computational efficiency becomes increasingly obvious.


Introduction
The geometric distribution of Global Navigation Satellite System (GNSS) tracking network has an important impact on the satellite orbit determination (OD), earth rotation parameter determination (ED) and geocentric movement monitoring (GM), among others. During the early GNSS system construction, only a few GNSS tracking stations can be can be obtained using the data from 8 stations within China's territory, but the OD precision can be improved markedly if an abroad station is added. The optimal location of such an abroad station is in Perth, Australia [Wen, Liu, Zhu et al. (2007)]. He found that BDS Geostationary Earth Orbit (GEO) satellite orbits, especially the along-track component, can be significantly improved by extending the tracking network in China along longitude direction, whereas Inclined Geosynchronous Satellite Orbit (IGSO) gain more improvement if the tracking network extends in latitude [He, Ge, Wang et al. (2007)]. Zhang analyzed the effects of ground tracking stations distribution on the precision of BDS satellite orbit determination, and proposed the optimal stations distribution strategy of BDS satellite orbit determination [Zhang, Dang, Cheng et al. (2016)]. In this contribution, the optimal design of GNSS ground tracking network is investigated and the determination precisions of satellite orbit and Earth Rotation Parameters (ERP) are considered simultaneously. And a method of rapidly determining stations optimal distribution locations based on map grid zooming and genetic algorithm is proposed. In the following, Section 2 presents a new criterion for the optimal stations distribution for OD and ED, named the minimum Orbit and ERP Dilution of Precision Factor (OEDOP) criterion; Section 3 presents a method for the rapid determination of the selected station locations based on the map grid zooming and genetic algorithm; Section 4 proposes an optimal or near-optimal distribution locations of 17 newly built BDS stations by using of the minimum OEDOP criterion and station optimal distribution location determination method; Section 5 summarizes the main conclusions and contributions of this paper.

Criterion for the optimal GNSS station distribution for orbit and ERP observation
The main goal of designing a GNSS global tracking network is to form a ground-based continuous observation network with the minimum number of stations for the precise determination of the satellite orbit, satellite clock, ERP, etc. However, in practical construction, the choice of sites has to consider a number of limitations, such as the observation environment, data communication conditions, electricity infrastructure, legality of site coordinates and so on [Hofmann-Wellenhof, Lichtenegger and Collins (2013)]. However, in this contribution, only two factors are considered: 1) The advantages of constructing a GNSS station on land relative to constructing it in the ocean and 2) The geometry structure of GNSS ground tracking stations for satellite orbit and ERP determination. The specific method is as follows. According to the principle of GNSS observation, the GNSS carrier phase observables are formulated as [Xu (2007)]: where  is the wavelength of carrier phase; Ф is the observed phase; tr is the signal reception time of the receiver i; te is the signal emission time of the satellite k; subscript i and superscript k are the number of receiver and satellite; tr and tk are the clock errors of the receiver and satellite at the time tr and te; c is the speed of light; N is the integer ambiguity; ion, trop, tide and rel are the ionospheric, tropospheric, tidal and relativistic effects, respectively; p is the remaining errors;  k i is the theoretical distance from the satellite k to the receiver i. From Eq. (1), it is known that the sum of left terms (the observational distance and all errors corrections) should be equal to the right term (the theoretical distance). However, the all errors cannot be completely corrected in the GNSS practical measurement. Thus, the following Eq. (2) is commonly used to represent the functional relationship between the GNSS observation distance and the theoretical distance. And a matrix expression is widely adopted to estimate the unknown parameters.
where L is the observation vector that includes the left terms of Eq. (1); A is the coefficient matrix that includes the all partial derivatives of Eq. (1) with respect to the unknown parameters and X is the unknown parameter vector; AX is the theoretical distance vector that includes the right term of Eq. (1); V is the residual vector that includes the differences between the observational distances and the theoretical distances; P is the symmetric and definite weight matrix that is used to assess the contribution or accuracy of each observational distance, which is commonly determined by the satellite elevation angle or signal-to-noise ratio of the each observational distance. However, in GNSS data processing, the priori values of unknown parameters are commonly used to obtain the partial derivatives of unknown parameters in the matrix A. Thus, the X includes the corrections of the priori values of the unknown parameters; the L includes the differences between the observational distances and the priori theoretical distances. And an iterative algorithm is used to obtain the optimal estimates of unknown parameters. It notes that the unknown parameter vector X just includes satellite orbit coordinates and ERP in this paper, and the other parameters are taken as the known values for the convenience of later discussion. The specific forms of the A, X, L and P are described as follows. Assuming that there are n GNSS stations and that each station tracks m satellites per epoch (the total number of epochs is j), A, X, L and P can be written as are the partial derivatives of the geometric distance ρ between the k-th satellite and the g-th station with respect to the x, y, z components of the satellite orbit coordinates; and are the partial derivatives of ρ with respect to Ɵx, Ɵy, Ɵu, where k={1…m}, g={1…n}. The formulas used to calculate these partial derivatives are as follows: Second, the inner coincidence of the i-th element of the estimated parameter in Eq.
(2) can be evaluated by where m0 is the so-called standard deviation (or sigma); pi is the i-th element of the precision vector, Qii is the i-th diagonal element of the matrix Q (the inversion of the normal matrix M); and where j, n and m are the consecutive natural numbers starting from 1. To directly describe the precision of a group of unknowns without knowledge of m0, a socalled dilution of precision factor (DOP) is usually used in GNSS measurements. The DOP is defined as the square root of the sum of some diagonal elements of the quadratic matrix Q, yielding the Position DOP (PDOP), Time DOP (TDOP), Geometric DOP (GDOP), etc. [Xu (2007)].
In this study, this definition is similarly used to describe the precision of OD and ED, where the corresponding DOP is termed the OEDOP. The calculation formula of the OEDOP is as follows: Generally, smaller the OEDOP values correspond to better the normal-equation condition of the estimated parameters. According to the above formulas, it can be readily noted that only the coordinates of the satellites and stations govern the OEDOP value. The satellites coordinates can be calculated by the precise or broadcast ephemeris and they cannot be changed since they have been defined in advance and controlled by the GNSS operation & control center. Therefore, the different OEDOP values will be obtained by using the different ground stations distribution scheme. Finally, comparing the OEDOP values for all possible station distribution schemes, the scheme with the lowest OEDOP value is selected as the optimal design scheme for satellite orbit and ERP determination.

A method for rapidly determining optimal station distribution locations
Using the above-proposed criterion of the minimum OEDOP value, the OEDOP values of all possible station distributions must be calculated and compared, which requires impractically high computation times. For example, if using a 20°×20° grid to cover the globe, 180 gird points are taken as potential station locations. However, 109 grid points are located in the ocean and are thus not suitable for GNSS stations. Therefore, the remaining 71 grid points located on land are taken as the final candidate locations. If we need to establish 17 new BDS stations, the total number of possible stations distribution schemes will be C71 17 =1.03×10 16 , where every scheme has a corresponding OEDOP value. It usually takes one hour to calculate 104 OEDOP values using a personal computer (e.g. a Dell MT3020). Therefore, it will take 1×10 14 hours (approximately 11.4 billion years) to complete all computations required to establish 17 new BDS stations among the 71 candidate locations using a 20°×20° grid. If a 1°×1° grid is used to improve the precision, the number of land-based grid points increases to 21616, which would require an inconceivably long computation time. Therefore, a method to rapidly select the distribution scheme with the minimum OEDOP criterion from all possible distribution schemes is urgently needed. To solve this problem, a new method is developed in this contribution. This method is inspired by the use of scale adjustment to find a certain location on a map. First, the scale of the map is decreased to quickly identify the region in which the location lies, and then the scale of this area is increased stepwise to find the specific location accurately. The merit of this method is that the location of interest is not compared with all the locations on the map. Based on this approach, a new method for the rapid determination of the optimal station distribution locations is developed. The specific implementation steps are as follows: (a) Determination of the general distribution First, a relatively sparse grid (e.g. 40°×40°) is used to cover the globe. Next, an initial optimal distribution of stations is obtained based on the given criterion (e.g. the minimum OEDOP). The goal of this step is to obtain a general configuration of optimal station distribution locations as quickly as possible.

(b) Densification of appointed grid
Second, a relatively dense grid (e.g. 20°×20°) is used to refine the grids nearest to the chosen grid point in Step (a). Next, all the grid points in the refined area are taken as the new candidate station locations. Fig. 1 shows the basic flows of Steps (a) and (b).

(c) Adjustment of chosen station locations
Third, the chosen optimal station locations are adjusted within the scope assigned in step (b) according to the minimum OEDOP criterion. It should be noted that the number of possible station distribution schemes is now Πi=1 n C25 1 (n is the number of stations), not Cn×25 n , because the general distribution has been determined in step (a). To reduce the computational time, a genetic algorithm is adopted, as usually used in the optimal design of GNSS networks [Saleh and Dare (2003); Saleh and Chelouah (2004)].

(d) Repetition of steps (b) and (c)
Step (c) yields new optimal station locations with lower OEDOP values than the locations found in Step (a). However, Steps (b) and (c) must be repeated if the grid density does not satisfy a predetermined criterion. These new optimal station locations are then taken as the new central points. The nearest grids are refined by a denser grid (e.g. 10°×10°), and the optimal station locations are adjusted again according to the minimum OEDOP criterion, as shown in Fig. 2. This cycle is repeated until the grid density reaches the required density (e.g. 1°×1°). In general, this method can markedly reduce the computation time and be easily automated and run using a personal computer. However, three questions must be answered to assess the effectiveness of the method: 1) Are the chosen optimal station locations obtained using the proposed method the same as those obtained using the exhaustive method? 2) How much can the calculation efficiency be improved using the proposed method relative to the exhaustive method? 3) Is it true that the smaller OEDOP values indicate higher precisions of the measured OD and ED? Lat.(

Experiments and analysis
To verify the above three problems, two sets of experiments were designed, and several sets of GNSS real data are processed. In the first set of experiments, the determination accuracy and computational efficiency of proposed method are tested. First, a 10°×10° grid is used to cover the globe, and the grid points located in the ocean are removed using the Matplotlib software and GSHHG high-resolution coastline data, which allow us to determine whether the points are on the mainland, an island or water [Janert (2011)]. Second, the exhaustive method and proposed method are adopted to determine the optimal locations for the construction of 1, 2 or 3 new BDS stations based on 13 existing stations using the minimum OEDOP criterion. Finally, the optimal station locations and computation time of the two methods are compared.
In the proposed method, the order of grid scales is 40°×40°, 20°×20°, and 10°×10°, and a genetic algorithm is used to reduce the computation time. In this genetic algorithm, the initial population is 10, and the probabilities of crossover and mutation are 0.6 and 0.4, respectively. These parameters are set based on the following rules: (1) 41 grid points can be generated if using a 40°×40° grid to cover the globe. However, there are 18 grid points are located in the ocean and 13 grid points are used to stand for 13 existing stations. Thus, 10 remaining grid points are taken as the alternative grid points, namely the initial population.
(2) The probability of crossover is calculated by the existing grid points divided by the total of the existing grid points and the alternative grid points, namely 13/(13+10)=0.6. (3) The probability of mutation is obtained by subtracting the probability of crossover from 1, namely 1-0.6=0.4. The BDS satellite coordinates for the two methods are obtained from the precision ephemeris of the GFZ GNSS analysis center, and Dell MT3020 computers are used in all experiments. Finally, only MEO and IGSO satellites are considered because the orbit determination precision of GEO satellites is far lower than that of MEO and IGSO satellites [Zhang, Zhang, Huang et al. (2015)]. Figs. 3(a), 3(b), 3(c) and 3(d) show the optimal station locations and computation time of the two methods.
(a) One optimal BDS station location (b) Two optimal BDS station locations (c) Three optimal BDS station locations (d) Computation time Figure 3: One, two and three optimal BDS station locations and computation time of exhaustive method and the proposed method From the above experimental results, several conclusions can be drawn. First, both methods yield the same optimal station locations in the cases of 1, 2 and 3 optimal BDS stations. However, the computation time of the proposed method is shorter than that of the exhaustive method. Furthermore, as the number of newly built stations increases, the proposed method offers an increasingly obvious improvement in computational efficiency relative to the exhaustive method. When the number of newly built stations reaches 3, the computational efficiency is 35 times better using the proposed method than using the exhaustive method. According to the BDS construction plan, 17 BDS abroad stations will be built in the next few years, which combined with the 13 existing stations, will form a global BDS tracking station network. However, this construction will be not completed within an acceptable time if the exhaustive method is used to determine the optimal locations of the 17 new BDS stations. Therefore, the proposed method is used instead. The order of grid scales is 60°×30°, 30°×15°, 10°×10°, 5°×5°, 2°×2°, and 1°×1°, and the initial population, crossover probability and mutation probability of the genetic algorithm are 10, 0.6 and 0.4, respectively. It takes approximately 18 hours to complete all computations using a DELL MT3020 computer. Fig. 4 shows the determination results in the initial 60°×30°, the intermediate 30°×15°, 10°×10°, 5°×5°, 2°×2°, and the final 1°×1° grid stages. From Fig. 4, it is known that the locations of 17 BDS abroad stations are refined step by step from the initial locations in Fig. 4(a) to the final locations in Fig. 4(f), which are obtained by using the proposed method that is described in Section 3. And it is noted that the adjustment of each station location is decreasing when the mesh size of the grid becomes smaller and smaller. This illustrates the optimal location of each station is obtained by the successive approximation method.
(e) 2°×2° (f) 1°× 1° Figure 4: Evolutional process of determining the optimal locations of 17 BDS stations using map grid zooming method and genetic algorithm Next, a second set of experiments is designed to verify that smaller OEDOP values indicate higher precisions of the OD and ED. These experiments include four experimental schemes and utilize real GPS observation data. All 30 stations are considered in each of the four experimental schemes, but the station distributions differ (see Fig. 5). Therefore, the OEDOP values of the four schemes are different. In Scheme 1, the stations are only distributed within Asia. In Scheme 2, the stations are distributed within Asia and Europe. In Schemes 3 and 4, the stations are distributed throughout the world, but the station locations in scheme 4 are from the results of Fig. 4(f). One week (from 2015.06.01 to 2015.06.07) of GPS data from these stations was gathered and processed, and one day of observational data was obtained as a solution session. High-precision GNSS data-processing software developed by Dr. Maorong Ge of the German Research Centre for Geosciences was used to process the experimental data [Li, Ge, Dai et al. (2015)]. The predicted values of IERS Bulletin A were used as the initial values of the ERP. The GPS broadcast ephemeris was used to generate the initial orbit. The coordinates of the station and their constraint values were obtained from the IGS SINEX file. The specific parameter configuration is listed in Tab. 1. Notes: The PC+LC indicates the ionosphere-free combination of pseudo-range and carrier phase; LSQ is the least square estimation; ISB is the inter-system bias and IFB is the inter-frequency bias. PCO and PCV are the absolution antenna phase center offset and phase central variation, respectively; SRP means solar radiation pressure.
The final ERP and GPS orbit productions of the IGS were taken as the reference values to evaluate the ERP and orbit solution precisions of each scheme. The differences between the IGS productions and the experimental results are shown in Fig. 6.  From Fig. 6, some conclusions can be drawn: 1) There are obvious systematic errors in scheme 1 when the observations only come from Asian stations, especially for the ORB and dXpole determination (e.g. Figs. 6(a) and 6(d)). The primary reason is that these unknown parameters are very sensitive to the distribution of stations, i.e. a minor change in the distribution of stations will lead to a great impact on the accuracies of ORB and dXpole parameters estimation; 2) the ERP and ORB determination precisions of Scheme 2 are better than those of Scheme 1 because the stations are more widely distributed. In particular, the systematic errors of Ypole are substantially removed; 3) the results of Schemes 3 and 4 are better than those of Schemes 1 and 2 because the stations have a global distribution, and the solution precisions of Scheme 4 are better than those of Scheme 3. Thus, it is verified that the station distributions in Fig. 4(f) are the optimal locations of 17 BDS abroad stations.
Tab. 2 lists out the statistical results from the different experimental schemes, and some conclusions can be obtained. 1) The accuracy of ORB determination from the scheme 4 is the best, which is improved by 69.1%, 38.5% and 12.6%, respectively compared with the results from the scheme 1, 2 and 3.
2) The accuracy of ERP determination from the scheme 4 is also the best, which is improved by 57.1%, 46.9% and 32.4%, respectively compared with the results from the scheme 1, 2 and 3. 3) And the results clear verify that smaller OEDOP values correspond to higher precisions of OD and ED.

Conclusions
During the early stage of GNSS system construction, only a few ground stations can be established due to organizational and financial constraints. Even if the number of stations is sufficient, only GNSS data from a few stations are used to reduce the computation time in some real-time applications. Therefore, optimizing the station distribution is a key issue, as it will allow the OD, ED and GM to be measured at a desired precision using a limited number of stations. In this contribution, the station optimal distribution problem is discussed in detail, using the optimal distribution of the 17 remaining BDS stations as an example. To judge the OD and ED precision of the different distribution strategies, an index value called OEDOP is proposed and tested. Then, to improve the computational efficiency of determining the specific optimal station locations, a fast method is developed based on the map grid zooming and genetic algorithm. The results show the computational efficiency of the proposed method is 35 times greater than that of the exhaustive method in the case of 3 newly built stations. And as the number of stations increases, the improvement in the computational efficiency becomes increasingly obvious. At present, constructing multi-GNSS stations is a global trend. Thus, it is important to consider which station distribution is most beneficial for the satellite orbit determination of multi-GNSS systems. Different GNSS systems have different satellite orbital inclinations and orbital periods, which have different sensitivities for ERP determination [Lutz, Meindl, Beutler et al. (2013); Rothacher and Weber (1999); Yang, Li, Xu et al. (2011)]. Therefore, which station distribution is optimum in terms of being most suitable for ERP determination using multi-GNSS observation data is also an important topic of investigation. It should be investigated in our future study.