Mathematical Problems in Engineering

Volume 2019, Article ID 7579185, 17 pages

https://doi.org/10.1155/2019/7579185

## Random Optimization Algorithm on GNSS Monitoring Stations Selection for Ultra-Rapid Orbit Determination and Real-Time Satellite Clock Offset Estimation

^{1}NASG Key Laboratory of Land Environment and Disaster Monitoring, China University of Mining and Technology, Xuzhou 221116, China^{2}School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China^{3}Satellite Positioning for Atmosphere, Climate and Environment (SPACE) Research Centre, School of Science, Mathematical and Geospatial Sciences, RMIT University, Melbourne, VIC 3001, Australia^{4}Chinese Academy of Surveying and Mapping, Beijing 100830, China

Correspondence should be addressed to Qianxin Wang; nc.ude.tmuc@xqw

Received 4 December 2018; Accepted 29 January 2019; Published 14 February 2019

Academic Editor: Antonio Elipe

Copyright © 2019 Xu Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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