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

Qianxin Wang , Chao Hu and Ya Mao

Abstract: 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.

Keywords: Global Navigation Satellite System (GNSS), optimal distribution of station network, map grid zooming, genetic algorithm.

1 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 available for organizational and financial reasons. Therefore, optimal design of the tracking network is a key issue, as it allows the determination of OD, ED and GM at the desired precision using a limited number of stations [Dvorkin and Karutin (2013)]. On the contrary, in the maturity stage, a number of ground tracking stations, GNSS satellites and signal frequencies can be available, which makes the burden on data processing increase significantly. In order to minimize the data processing burden without significantly sacrificing the precision of OD and ED, the optimal minimum stations need to be identified and selected, particularly in the area of real-time and near real-time applications [Wang, Zhang, Wu et al. (2017)]. No matter in the initial stage or in the maturity stage, the optimal design and optimal utilization of GNSS stations are always very important for the parameters estimation precision and the computational efficiency.Therefore, many investigations have been performed to tackle this problem [Wang, Dang and Xu (2013); Zhang, Zhang, Huang et al. (2015); Hu, Wang, Wang et al. (2017)].

During the early development of GPS, the optimization of GPS global tracking network design was studied widely for various purposes [Delikaraoglou (1985); Bernd and Wolfgang (1988); Stephen and Willy (1989)]. These studies have found that longer distances between monitoring stations yield better precisions of satellite OD. Furthermore,their results indicate that a minimum of 15 globally distributed stations is required for centimeter-level orbit determination, beyond which the improvement provided by additional stations goes as the square root of the number of stations [Wang (1997)]. With regard to GLONASS, GALILEO, the optimal designs of their respective global telemetry,tracking and command (TTC) networks have also been investigated. Dvorkin studied the optimization distribution of GLONASS global tracking stations network. The results show the precision of satellite orbits and clocks (OC) can be improved significantly by adding 11 abroad stations to the 10 existing stations, which provide 4-fold coverage.However, increasing the number of stations by almost 30% when establishing a network with 5-fold coverage, it increases the accuracy of OC by 0.5% [Dvorkin and Karutin(2013)]. Liu discussed the optimal design of GALILEO tracking network and verified that the accuracy of satellite OD can be improved by about 80% when adding 9 stations to the 12 Cooperative Network for GIOVE Observation (CONGO) stations. The new monitoring stations are most effective when located in North America, Asia, the Arctic and Antarctica [Liu, Wu, Cai et al. (2003)].

The BeiDou Navigation Satellite System (BDS) is a global navigation satellite system,which has been independently constructed and operating by China. The deployment of BDS is divided into two phases: the regional system and the global system. By 2012, the regional system has been constructed, and the global system has been developed since 2015. According to the BDS construction plan, by 2020, 30 BDS satellites will be launched as a constellation, and about 30 ground stations (8 domestic and 22 abroad) will be built to form a global BDS tracking network [CNSO (2012)]. As of 1st January 2015,14 BDS satellites have been launched successfully, and 13 BDS tracking stations have been established, which can meet the service requirements for China and its surrounding areas [Guo, Li, Zhang et al. (2017)]. However, in order to provide precise position,navigation and timing (PNT) service to global users, the optimal design and construction of the BDS global ground tracking network is very important. And some investigations of this problem have been performed. Wen demonstrated that sub-meter BDS satellite OD 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.

2 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; tris the signal reception time of the receiver i; teis the signal emission time of the satellite k; subscript i and superscript k are the number of receiver and satellite; δtrand δtkare the clock errors of the receiver and satellite at the time trand te; c is the speed of light; N is the integer ambiguity; δion, δtrop, δtideand δrelare the ionospheric, tropospheric, tidal and relativistic effects, respectively; εpis the remaining errors; ρkiis 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

Here, Xsatand Xerpare the unknown parameter vectors of orbit and ERP, where Xsat=[x, y,z]T, Xerp=[?x, ?y, ?u]T, and x, y and z denote the coordinates of the satellite orbit on the X,Y and Z components. Furthermore, ?x, ?yand ?udenote the pole motion on the X, Y components and the value of ut1-utc, respectively. Asatand Aerpare the coefficient matrices of Xsatand Xerp, respectively. Moreover, matrices Asatand Aerpat the i-th epoch(i={1…j}) can be written as

Second, the inner coincidence of the i-th element of the estimated parameter in Eq. (2)can be evaluated by

where m0is the so-called standard deviation (or sigma); piis the i-th element of the precision vector, Qiiis 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.

3 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 C7117=1.03×1016, 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×1014hours (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=1nC251(n is the number of stations), not Cn×25n, 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)].

Figure 1: Diagram of Steps (a) and (b) for rapidly determining optimal station distribution locations

(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°).

Figure 2: Diagram of steps (c) and (d) for rapidly determining optimal station distribution locations

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?

4 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.

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.

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).

Figure 5: Different stations distribution locations for Schemes 1-4

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.

Troposphere Dry and wet GMF mapping functions; ZTDs estimated for each station at intervals of two hours Precession and nutation IAU 2010 precession and IAU 2010 nutation model Geopotential EGM2008 12×12 Solid Earth tides, ocean tides and solid Earth pole tides IERS Conventions 2010 N body JPL DE405 ephemeris Orbital parameters Six orbit elements and five ECOM SRP parameters were estimated: constants in the D, Y and X directions; periodic terms in the X direction Pseudo-stochastic orbit parameters Every 12 h; constrained to 1×10-6 m/s in the radial direction, 1×10-5 m/s in the along direction, and 1×10-8 m/s in the cross-track direction Ambiguity Real constant value for each ambiguity arc

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. The average values of OEDOP in this week and the practical solution precisions of the four schemes are provided in Tab. 2.

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.

5 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.

Acknowledgement:This work was supported by “The National Natural Science Foundation of China (No. 41404033)”, “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 Jiangsu Dual Creative Teams Program Project Awarded in 2017” and thanks for the data from IGS and iGMAS.