A Navigation Satellites Selection Method Based on ACO With Polarized Feedback

Selecting the optimal satellite subset for positioning from all satellites in view can not only achieve positioning accuracy but also reduce the computational burden. In this article, a navigation satellites selection method is proposed based on ant colony optimization with the improvement of polarized feedback (ACO-PF). Firstly, the satellite selection problem is described as a combinatorial optimization problem, and the noise weighted geometric dilution of precision (NWGDOP) is defined as the criterion for satellite selection. Then the ant colony optimization (ACO) is incorporated to solve the problem, and a polarized feedback mechanism is presented to improve the convergence speed of algorithm. Meanwhile, a perturbation operator is designed to improve the global searching ability of the algorithm. The numerical experimental results show that ACO-PF can select the superior satellites combination which provides high-precision positioning. And its convergence outperforms the related algorithms by up to 50%. Besides, the achieved NWGDOP of ACO-PF is usually 0.065 smaller than ACO. Therefore, the ACO-PF method can be considered as a promising candidate for satellite selecting in navigation applications.


I. INTRODUCTION
The Global Navigation Satellite System (GNSS) is the system that provides positioning, navigation and timing (PNT) services for the global region [1]. GNSS includes the well-developed Global Positioning System (GPS) of USA, GLONASS of Russia, and BeiDou Navigation Satellite System (BDS) of China, Galileo satellite navigation system of European in active construction. With the gradual expansion of their coverage, the applications and services relying on the GNSS are also increasing [2]. The GNSS plays an increasingly important role in our work and daily life. How to effectively use GNSS to serve us is a necessary and meaningful task. The selection of navigation satellites combination is an extremely important part of it. At present, in addition to The associate editor coordinating the review of this manuscript and approving it for publication was Massimo Cafaro . navigation and positioning in single system, multiple systems integrated navigation and positioning is also developing.
As for multi-system observation space, the number of visible satellites multiplies. With the deployment of Galileo and Beidou systems, and the modernization of GLONASS, there will be four GNSS constellations providing totally more than 40 visible satellites [3]. Because the unknown number of synchronization time deviation among the satellite systems increases, more satellites need to be selected for positioning [4]. In practical applications, we need to use the navigation messages of satellites to compute the position of the receiver. But if all observable satellites are used to compute the solution, the receiver will encounter extremely burdensome pressure due to large computational complexity [5]. In order to guarantee positioning precision while reducing computation load, it is necessary to select the right satellites subset from observable ones for solving problem in the various applications mentioned above. The problem of satellites selection has been attracted many academicians to study. The key to solve this problem is to consider the positioning precision and real-time performance comprehensively. In this way we can design reasonable and effective satellites selection algorithm.
Previous researches have found that a performance metric of selected satellites subset that could reflect the 3D geometric configuration of satellites, which is represented as the geometric dilution of precision (GDOP) [6]. It generally think that the smaller the GDOP is, the more accurate the navigation fix is [7]. Based on the above findings, the minimum GDOP value (MGV) method, the largest tetrahedron volume method and the maximum orthogonal projection were found and widely used. As the traditional methods [8], the computational load increases significantly with the number of observable satellites increases. Besides, using GDOP as criterion for satellite selecting only considers the satellites geometry. In fact, the positioning precision is affected by many other factors, such as the error of ephemeris from the space vehicle, the multi-path effect from the signal transmission path and the SNR from GNSS receiver [9]. The latter two methods are mainly used for exactly four-satellite subset selecting. So these methods are usually not applied to engineering projects directly. In recent years, a lot of researchers proposed some machine learning (ML) methods to approximate and classify GDOP value to solve the satellites selecting problem [10]- [13]. In addition to the shortcomings of these methods themselves such as easy to trap into local optimum, it is also difficult to implement large-scale training samples. Besides, some satellites selecting methods using neural network (NN) usually give the approximate value of GDOP and their performance also depend on training process as different training patterns affect the performance differently [12], [13]. Different from the methods using brute force selecting procedure, we consider satellites selection as an objective combinatorial optimization problem and adopt an intelligent optimization method ACO-PF to search the optimal solution. The improvements for the method can greatly improve the algorithm performance. At the same time, the criterion for selecting satellites in ACO-PF is noise weighted GDOP (NWGDOP) which considers the multi-path errors caused by satellite signal strength and noise. The improved satellite selecting criterion can further improve the performance of selected satellite combination.
The rest of the paper is organized as follows. The related research work in recent years was analyzed and summarized in Section II. The problem of navigation satellites selecting is specifically represented with the steps to obtain NWGDOP value in Section III. The ACO-PF method is presented in details and applied to solving the satellites selecting problem in Section IV. Extensive experimental results are provided in Section V before conclusions are drawn in Section VI.

II. RELATED WORK
As mentioned above, GDOP is one of the main performances to judge the quality of selected satellites subset. The original satellites selection method is minimizing the GDOP value of satellite subset, which is called the minimum GDOP value (MGV) method [14]. By this method, all possible satellites subsets are selected from all observable satellites, then select the one with the minimum GDOP value as the optimal subset [15]. This is a simple permutation and combination calculation. As the number of satellites increases, it will form a combinatorial explosion problem and the computational complexity will also increase to an unacceptable level. So the MGV method is just a theoretical method of satellite selection which is suitable for small number of satellites in single constellation. To relieve the calculation burden of this original method, there are some quasi-optimal GDOP methods proposed.
A. Peng et al. proposed a recursive satellite selection method to expand the selected satellites subset sequentially [16]. However, the polynomial time complexity of this method is worthy of improving. F. C. Meng combined fast satellite selection algorithm and receiver autonomous integrity monitoring (RAIM)/ fault detection and exclusion (FDE) by weighted matrix estimation in weighted least squares (WLS) [17]. But this method is also the estimation method of GDOP value. M. M. Wei et al. proposed a new algorithm that makes use of two satellites' comparability and the distribution characteristic of satellites with high or low elevation to select satellites [18]. D. Roongpiboonsopit et al. proposed a multi-constellations satellite selection algorithm (MCSSA). The MCSSA considers the geometric arrangements of observable satellites and selects the ones that spread evenly over the observed sky view [19]. M. Zhang et al. put forward an idea that selecting a subset of all satellites in view whose geometry is the most similar to the optimal geometry as the optimal satellite combination [20]. However, these methods usually select the sub-optimal satellite subset merely. Besides, the positioning precision of the navigation fix is usually not satisfying. In [21], H. Liang et al. ensured the positioning accuracy of the satellite from a new point of view which is by reducing the low elevation multi-path interference signal. However, the determination of the optimal mask angle is based on the DOP value, which also needs to be deduced by the author.
Considering the disadvantages of the MGV method with brute force, some academicians put forward geometrically optimized satellites selection algorithms including maximizing the volume of sagittal tetrahedron of the satellites subset [22]- [25]. This method selects the optimal satellites subset based on the relations between GDOP value and geometric form of selected satellites. Similarly, another method named the maximum orthogonal projection is proposed [26]. However, the two methods based on the geometric form of selected satellites are only suitable for selecting four satellites from all observable satellites, which can only meet the needs of VOLUME 8, 2020 satellite selection in single constellation. Through the two methods, as the number of satellites increases, the location precision will decrease and the amount of calculation will increase. Many scholars have begun to study alternative methods to GDOP computation.
In recent years, ML-based methods such as NNs, genetic algorithm (GA), and support vector machine (SVM) algorithm are used to solve satellites selection problems. Almost all satellite selection optimization methods based on ML algorithm are mainly divided into two categories, that is, approximation and classification. A group of methods based on approximate calculation are proposed in order to simplify the calculation of GDOP, such as NNs [27], SVM [28]. However, the GDOP approximation methods just treat satellite selection as a regression problem of GDOP calculation, and the optimal subset is also need to be selected from all satellite subsets by the original brute force means. Besides, another group of machine learning (ML) based methods focus on classifying GDOP value for selecting one of the acceptable subsets for navigation use. Dan Simon and Hossny El-Sherief have used the neural network (NN) approach create a network which classifies satellite subsets according to GDOP as early as 1995 [12]. N. Zarei used an artificial neural network (NN) that can generate GDOP without inverting a matrix [29]. It is based on the learning properties of artificial neurons, but as such gives an approximate rather than an exact answer. Since the performance of some of these methods also depends on the training time and the size of the training data, D. J. Jwo et al. presented the neural network (NN)-based navigation satellites subset selection to reduce the training time of back propagation neural network (BPNN) [30]. A. Hamed et al. proposed an improved BP algorithms, including resilient BP that has greater precision and calculation time [31]. Even adopting classification thought, these methods can only approximately classify the GDOP value of each satellite subset into a predefined range [32].
Evolutionary algorithms are also been researched for using satellite selection. In 2010, M. R. Mosavi proposed using GA to achieve satellite selection, and reduced the number of calculations of GDOP through the rapid optimization of GA [10]. Since GA result lacks hill climbing ability and is easy to trap into local optimum in the later iterations [33], in 2012, M. R. Mosavi et al. discussed the GPS satellites classification algorithms based on evolutionary algorithms (EAs) [34]. Using EAs in the classification process avoids calculating inverse matrix so that it can reduce the calculation burden of the computation process. However, since EAs do not have a clear genotype-phenotype distinction, they are not reliable in many applications such as GPS GDOP classification [13]. Besides, only the GPS system is involved in [33], which has limitations. The researchers have made a lot of improvements to the GA method, which improves the precision of the algorithm [11], [35]- [37]. However, all these research is also focus on GDOP classification by GA-based methods. Besides, the excessive number of parameters need to be adjusted in the calculation process of GA.
In addition, most of the above methods are based on GDOP value as a satellite selection criterion. Since the GDOP value represents the geometric relationship between the observation point and the satellites for navigation use, it doesn't consider the error statistics property which causing by received signal of satellite and interference noise in space.
In order to avoid the shortcomings of the above existing methods, this article proposes a navigation satellites selecting method based on ant colony optimization with polarized feedback (ACO-PF) to effectively improve the timeliness of navigation satellites selection and the positioning precision of the selected navigation satellites combination in single constellation or multi-constellations. Since the GDOP value represents the geometry of the satellites for navigation use, we considers the signal-to-noise ratio (SNR) of the satellite on the basis of GDOP, thus forming the NWGDOP criterion.

III. PROBLEM FORMULATION
In general, selecting the optimal satellites subset can be classified as combinatorial optimization problems; selecting N objects from a choice of M allows for C N M potential subsets that makes the objective function value as small as possible. Based on the above principle, we propose the following problem model.
In this article, the observable satellite refers to the satellite which can receive its signal at the observation point, and whose elevation angle is higher than the masking angle. Considering that when the elevation angle is less than 10 degrees, it is greatly affected by the atmospheric delay error [38], so the masking angle is set to 10 degrees in this article. Set S = {1, 2, . . . , M } as the index set of M observable satellites, s represents the index of each satellite in set S, so s ∈ S. Now we have the geocentric coordinates of M observable satellites Q = [X, Y, Z] and corresponding SNRs R. And define p = [x, y, z] as the geocentric coordinate of the observation point. The coordinate of each observable satellite is relative to the observation point. One satellite subset T which includes N satellites is selected from S. The satellite optimized selection problem is to search the optimal T from S using the coordinates and SNR informations.
where J is the minimum number of satellites to ensure the precision of satellite positioning, usually is 4 in single constellation and 5 in multi-constellation. And G function represents the objective function for given Q,R, p and T. Since Q and p are known, the G function value can be computed by T shown as following definition.  Satellites geometry distribution and positioning error: (a) more concentrated satellites distribution geometry with higher GDOP; (b) more evenly satellites distribution geometry with lower GDOP. The solid circles which represent the true distance from satellite to user intersect at receiver's position. The dotted circle parallel to the solid circle represents the erroneous measurements obtained by the receiver. The red field is the uncertainty region.
In GNSS systems, the accuracy of satellite positioning depends on many factors including the number of visible satellites, the geometry of the satellite combination and the condition of the visible signals [39]. The standard deviation of navigation error σ A can be represent as the product of GDOP and user ranging error: where σ UERE is the user equivalent range error (UERE).
From (2), we can obtain that the GDOP value and UERE determine the positioning accuracy. GDOP can be seen as the magnification of the total positioning error to the UERE [40]. Since the distribution of UERE resembles the Gaussian distribution and is equivalent for each satellite, GDOP can reflect the positioning accuracy of the satellite combination to a certain extent. GDOP represents the geometry of the satellite combination as seen from the receiver. The visual description of GDOP is shown in Fig. 1. In Fig. 1, the circles at the bottom of the figure denote the relation between positioning error and satellites geometry. The uncertainty region is larger as the satellites spatially scattered evenly in Fig. 1(a). The concentratedly distributed satellites have larger uncertain region as Fig. 1(b) shows. That is, the positioning error of a better geometry is smaller.
As the number of visible satellites is N , the observation matrix of N satellites is: where e i = [e xi , e yi , e zi ] is the unit vector of position of the receiver to the position of the ith satellite. And the definition of GDOP is: where tr(·) is the symbol of matrix trace computing.
From (3) and (4), it can be seen that the GDOP is solely dependent on the geometrical distribution of users and visible satellites in space. Taking GDOP as the satellite selection criterion, only the geometry distribution is considered. In fact, the satellite positioning precise is influenced by not only the satellites geometry but also the other factors such as the error of ephemeris, the multi-path effect from the signal transmission path and the SNR from GNSS receiver [41]. The analysis on these factors is complex and not clear enough. For example, considering only the geometry of the satellite combination and receiver may miss some satellites with valuable signals. Another example is that the satellite with lower elevation can effectively enhance the geometry of satellites combination, but it is easy to be interfered by multipath effect. This article supposes that the satellite positioning precise is influenced by not only the satellites geometry but also the elevation angle and SNR of satellite.
On the one hand, we proposed the masking angle to remove the satellites low elevation angle in advance. It is known that the signal from a satellite with a lower elevation angle is much more affected by multi-path than that from a higher elevation angle [42]. The multi-path error is relatively small when satellite masking angle is 10 degree thus the satellite masking angle is set to 10 degrees [38]. The masking angle is set to 10 degrees, that is, the satellites whose elevation were no less than 10 degrees are used as effective satellites. On the other hand, considering the impact of satellite signals and noise on positioning accuracy, we add noise weighting matrix W to the definition of GDOP. W is a diagonal matrix and defined as a weighted matrix: where w i is the SNR of ith satellite in set S. The NWGDOP is defined as follows: SNR refers to the ratio of received carrier signal strength to noise strength, which can better reflect the quality of satellite signal. The higher the SNR value is, the better the carrier signal quality is. Based the definition of NWGDOP, we can conclude that the characteristics of NWGDOP. Firstly, the satellite with low SNR can be excluded directly to avoid affecting the positioning accuracy. Then the satellite with good signal quality but can not form good geometry with other satellites is considered. That is to say, NWGDOP criterion for satellite selection achieves a balance between satellite subset geometry and satellite signal quality. This article aims to minimize that NWGDOP value of the selected satellite subset. NWGDOP of one satellite combination is computed with the detailed information of selected satellites and observable satellites according to Section III-A and III-B. Since W is a diagonal matrix, the increased computational head of NWG-DOP compared to GDOP can be negligible.
On the other hand, the ACO algorithm is applied in more fields now and it will be applied frequently in combinatorial optimization problems in the future [43], [44]. K. F. Doerner et al. have researched the application of ACO in portfolio optimization problem [45]. X. S. Lai et al. applied it to fast matrix multiplication [46]. This article proposes an optimized ACO method, ant colony optimization algorithm with polarized feedback (ACO-PF), to select the optimal or nearoptimal subset from the set of all observable satellites. Back to the definition of GDOP in Definition 1, its detailed calculation procedure is introduced.

A. STEPS TO ESTABLISH NWGDOP FUNCTION
Set the observation point as the origin, the long axis of the earth ellipsoid as the X axis, the short half axis of the earth ellipsoid as the Y axis, and the normal of the earth ellipsoid as the Z axis, thus the carrier coordinate system is constructed [47].
where T is the transformation matrix between the carrier coordinate system and the geocentric coordinate system: where B p , L p are the latitude and longitude of the observation point that are known. After the above steps we got the carrier coordinates of each observable satellite. The relationship between the coordinates x c s y c s z c s T and the azimuth angle A s and elevation angle E s of the satellite is as follows: The azimuth angle A s and the elevation angle E s of the satellite s are obtained as following: Let the selected satellite subset be T = {S 1 , S 2 , · · · , S N }, then the azimuth and elevation angles of the navigation satellites are A S 1 , A S 2 , · · · , A S N and E S 1 , E S 2 , · · · , E S N respectively; the observation matrix H of the navigation satellite combination in single constellation is obtained: Considering the impact of noise, we construct the noise weighting coefficient of satellite subset be T as follows: We substitute the SNRs R of the satellite subset into (5) to get the noise weighting coefficient. Then we set NWGDOP function as fitness function of this paper's algorithm, which shows as follows: where the f (T) represents the NWGDOP value of satellite subset T, tr(·) function denotes the trace of the matrix; according to the research, the smaller the fitness function value f (T), the better the performance of the satellite combination T.
In the case of multi-constellation integrated navigation, the calculation of NWGDOP value needs to consider the time error of different systems so that the state matrix H in (12) should be revised accordingly. Taking GPS/BeiDou dual system integrated navigation as an example, H is represented as [48]: where A S i , E S i (i = 1, 2, · · · , N ) are the azimuth and elevation angle of satellite S i respectively, N 1 is the number of available satellites in GPS and the number of Beidou's satellites is (N −N 1). In the same way, the observation matrix of three-system or even four-system integrated navigation can be derived. Substituting (14), as shown at the bottom of the next page, into (13) can get the NWGDOP function in the case of multiple constellations.

B. THE SEARCHING SPACE OF SATELLITES SELECTING
In Fig. 2, the schematic diagram of 2-dimensional space with M satellites in each column is shown. There are N columns in searching space. The S i represents the ith satellite in set S. And Fig. 3 is the formalistic description of searching space defined as follows.
Definition 2 (satellite Searching Space P M ×N ): As illustrated in Fig. 3, the formalistic description of searching space including N columns is constructed, and each column contains M observable satellites; let O ip denote the ith observable satellite node in the pth column with i = 1, 2, . . . , M ; and denote O jq the jth observable satellite in the qth column node with p = 1, 2, . . . , N ; Let ants search their path in P M ×N and we called it satellite searching space.
Each ant selects a satellite in each column of the searching space P M ×N to form a navigation satellite combination with N satellites. Until all K ants complete the search, there will be K combinations. This search process is pushed to the  Through the above analysis, we can intuitively see that the navigation satellites selection problem is essentially an independent combinatorial optimization problem, which lays the foundation for us to solve it with ACO algorithm.

IV. PROPOSED ACO-PF ALGORITHM A. OVERVIEW OF THE ACO ALGORITHM
The ACO algorithm is a new type of simulated evolutionary algorithm firstly proposed by the Italian scholar Marco Dorigo et al [49], [50]. It is based on the study of the collective behavior of real ant colonies in nature. When ants search for food, they usually deposit a special chemical on the path they travel through. Other ants can perceive the pheromone concentration and guide their movement direction, so that they tend to move in the direction of high pheromone con-centration [51]. The ACO algorithm has been proved to be effective at finding optimal or near-optimal solution for some combinatorial explosive problems such as traveling salesman problem (TSP), allocation problem and job scheduling. Preliminary research has shown that it has innate advantages such as parallelization, positive feedback, and robustness in solving such complex combinatorial optimization problems [52], [53]. The navigation satellite selection problem is to find the optimal navigation satellite combination with the smallest GDOP value among all observable satellites (i.e., satellite searching space P M ×N ), which belongs to a class of discrete combination optimization problems. So this problem can be solved by the ant colony optimization algorithm. On this foundation, we propose two improvements to the basic ant colony optimization algorithm to solve the problem of navigation satellites selection, which presented as follows.

B. THE SPECIFIC STEPS OF ACO-PF FOR SELECTING THE OPTIMAL NAVIGATION SATELLITE COMBINATION
It can be seen from Fig. 3 in Sec. III-B that the satellite searching space P M ×N contains N columns with M satellites per column. Each ant selects a satellite in one column until all the columns are reached so that each ant obtains N non-repetitive satellites which format its moving path. After completing one iteration, K paths are formed by K ants if there are K ants in the ant colony. The algorithm repeats the above steps about searching paths until the maximum number of iterations is reached.
Denote l as the serial number of each iteration, let l max be the maximum number of iterations, and initialize l = 1. Let K be the total number of ants in the ant colony. The size of K is determined by experimental methods and experience, but the number of ants cannot be more than the total number of satellite nodes in the searching space. Therefore, K M · N . Define k as the serial number of ant and initialize k = 1. The set t k includes the satellite nodes that the kth ant has passed. Relatively, the set a k includes the satellite nodes that the kth ant never passed, that is a k = S − t k where S is the index set of M observable satellites defined in Sec. III. Set d k as the column corresponding to the initial observable satellite node selected by the kth ant.

1) ALGORITHM DESIGN BASED ON ANT COLONY OPTIMIZATION
According to the basic principle of the ACO algorithm, an ant can leave a kind of substance called a pheromone on the VOLUME 8, 2020 path, and the ant can perceive the information and move toward a high concentration of the pheromone. The ants are connected in this way, so the group's foraging behavior shows a positive feedback [54]. When ant k is in node O ip and has so far selected the partial satellites t k in the previous process, the probability of going to node O jq is given by: In above equations, τ ip,jq is the pheromone added to the edge (O ip , O jq ) on the current path, η ip,jq is a static greedy measure of the ''goodness'' of edge (O ip , O jq ) called heuristic information. Subscript i, j of the variables correspond to the satellite in each column in the searching space P M ×N and subscript p, q correspond to columns. O ju is a node not yet visited by the ant k. The parameters α and β control the relative importance of the pheromone τ ip,jq versus the heuristic information η ip,jq .
The ant is given priority to select the navigation satellite shown in (15) and (16) that can make the geometry of the navigation satellite combination as large as possible, so that the GDOP value is as small as possible [55]. That we can get the smallest fitness function value [56]. Therefore, the heuristic information η ip,jq is defined as the Euclidean distance between node O ip and node O jq : The pheromone τ ip,jq (t) on the connecting path between node O ip and node O jq is updated as (18) shows: where ρ ∈ (0, 1) represents a pheromone residual coefficient. τ k ip,jq represents an increment of pheromone on the path between the node O ip and O jq by ant k: where T k is the navigation satellite subset selected by ant k; Q is a positive parameter. The adjustment of the pheromone here uses the overall information, and the performance is better when solving the problem. Ant k finds the path with the probability shown in (15).

2) TWO IMPROVEMENTS FOR ANT COLONY OPTIMIZATION BASED ALGORITHM
The ACO algorithm is usually used to solve the difficult combinatorial optimal problems although it has the defects of slow convergence speed and is easy to fall into local minimum [57]. ACO-PF proposed in this article is an improved algorithm which contains two improvements added on ACO algorithm. One is the polarized feedback mechanism added on the corresponding variable. Another is the perturbation operator designed to affect one of the key parameters in the solving equation. The two improvements have guaranteed that the proposed algorithm performed superiorly, which highlighted in the Fig. 4. On the one hand, in order to improve the convergence speed of the algorithm, the calculation formula of pheromone increment is improved as (20), in which the polarized feedback mechanism is constructed: where T k denotes the navigation satellite subset selected by ant k; B is an empirical parameter, which is generally no more than 60% of the number of convergence iterations. The polarized feedback mechanism is introduced in the first B cycle of the algorithm and defined below.

Definition 3 (polarized Feedback Coefficient): ξ PFC is a polarized feedback coefficient to adjust τ k ip,jq according to the fitness function value. It is presented as follows:
where f (T k ) is the GDOP value of the navigation satellite subset T k selected by ant k; r is the GDOP reference value of navigation satellite subset T k , which is usually 6 according to experience. It can be seen from (21) that if the GDOP value of T k is smaller than r, the quality of the solution is indicated not bad, then the feedback coefficient ξ PFC is +1 so that the pheromone increment will get the corresponding ''positive'' increment. Otherwise, it means that the solution is not available, then the feedback coefficient ξ PFC is -1 so that the pheromone increment will get corresponding the ''negative'' increment. The above positive and negative polarized feedback mechanism can promote the superiors surviving, which is conducive to the rapid convergence in the initial stage of the algorithm (pre-B cycle). It can be seen from (21) that the polarized feedback mechanism is cancelled after the B loop, and the algorithm enters into the steady state convergence process. On the other hand, in order to improve the global searching ability and searching speed of the algorithm, a ''disturbing impact'' is introduced to make the algorithm jump out of the local optimal solution.
Definition 4 (Perturbation Operator): Based on the probability formula of ACO shown as (15) and (16), an influence factor d called perturbation operator is added to the parameter β. At the end of each iteration, d is adjusted according to the obtained optimal solution. Let T l be the optimal satellite combination of the lth iteration. In the lth iteration, K ants get K satellites subsets and T o l is the satellite subset with the minimum GDOP value among the K subsets. Then the perturbation operator d is updated to where d max represents the maximum value of d. Therefore, the probability P k ip,jq (t) of ant k from node O ip to node O jq is: As mentioned after (15) and (16), τ ip,jq represents the amount of pheromone on the edge (O ip , O jq ) in (23) and (24) similarly. The parameter α represents the relative importance of the pheromone quantity τ (ip,jq) ; η ip,jq represents the heuristic information. Both of the above formulas are satisfied q ∈ a k . Clearly, the parameter β represents the relative importance of the heuristic information η ip,jq . perturbation operator d is added to adjust parameter β to achieve the purpose of adjusting probability P k ip,jq adaptively. Combined with (22), (23) and (24), we can see that in the initial period of the algorithm operation, the optimal solution obtained is still in evolution. At this time, d is 1, which does not affect the probability P k ip,jq . To be further evolved, when the optimal solution obtained by the algorithm is not significantly improved, the perturbation operator begins to play a role. With the updating of perturbation operator d, the role of heuristic information η ip,jq in P k ip,jq is significantly enhanced, while the role of pheromone τ ip,jq is reduced. Thus, the ants tend to explore new solutions, which makes it easier for algorithm to jump out of local optimal solutions. At the same time, in order to ensure the convergence of the algorithm, the perturbation operator is restricted by d max . Once the solution jumps out of the local optimum, that is the optimal solution starts to evolve again, then d = 0, the perturbation disappears. Under the influence of perturbation operator, the solution jumps out as soon as the algorithm occurs in a local minimum, thereby continuing the evolution to the optimal solution.
In every iteration, ant k seeks forward and backward from the initial satellite nodes respectively, avoiding the algorithm falling into local optimum and promoting the evolution of the global optimal solution.

Algorithm 1 Algorithm Steps of ACO-PF Method
Init ← Let Pos k be the initial position of ant k; Denote C k as the number of column of Pos k ; for each k ∈ [1, K ] do p = C k ; while p < M do /*Select the node of latter one-column*/ O jq = arg(max(P k ip,jq ))(j ∈ a k ); add the number i into the set t k ; add the node O ip to the integer vector T k ; q = p; p + +; p = C k ; while p > 1 do /*Select the node of former one-column*/ O jq = arg(max(P k ip,jq ))(j ∈ a k ); add the number i into the set t k ; add the node O ip to the integer vector T k ; q = p; p − −; Update the pheromone τ k ip,jq according to Algorithm 2; Get local optimal solution T o l in lth iteration; Get global optimal solution from T o Compute the fitness function value set f (T 1 ), f (T 2 ), · · · , f (T K ) of the K navigation satellite combinations of the lth iteration completed by K ants using (7) to (12), and the local optimal navigation satellite of the lth iteration combination T o l is obtained by (25): where k ∈ {1, 2, · · · , K }; minf (T k ) denotes the minimum value of the fitness function value set of the K navigation satellite subsets. When the number of iterations reaches the maximum l max , we can get the global optimal navigation satellite combination T o which is smallest in all local optimal solutions. VOLUME 8, 2020 Based on the discussion, we present the complete navigation satellites selecting method in Algorithm 1 and 2. The main flow of ACO-PF method in this article is presented in Algorithm 1. And Algorithm 2 focuses the two improvements in ACO-PF method versus ACO method.
The main process of ACO-PF is show as Fig. 4 and its corresponding pseudo code description is present as Algorithm 1. Two improvement measures added to the basic ant colony algorithm, which are presented in Algorithm 2. One is perturbation mechanism applied to the desired heuristic factor, in order to prevent the algorithm from falling into local optimum. The other is the introduced polarized feedback mechanism, which aims to speed up the convergence of the algorithm.

C. COMPLEXITY ANALYSIS
From the algorithmic steps of ACO-PF in Algorithm 1, it can be seen that the worst case complexity of ACO-PF is O(NC · m · n 2 ). NC is the maximum number of iterations of the algorithm. m is the number of ants. n is the number of observable satellites.
The time complexity of the algorithm increases as the size of the problem increases and decreases as the number of ants decreases. The increased complexity due to the proposed improvements for ACO-PF method can be neglected, because it is insignificant compared to the main steps.

V. EXPERIMENTAL RESULTS AND ANALYSIS
As Fig. 5 shows, we use the high-precision GNSS receiver and measuring antenna with GPS/Beidou system, which can output the satellite navigation messages, to conduct the satellites selecting experiment. The data output frequency of the GNSS receiver is up to 20Hz, and the data is sent to the computer through RS232 port. The related methods are applied in the experiment using VC++ software.
As shown in Fig. 6, the GNSS receiver could receive signals from 21 satellites including 10 GPS satellites and 11 Beidou satellites. Considering the clock error of GPS system and Beidou system, the number of satellite selected is at least 5. And in the dual-system, properly adding the number of selected satellites is beneficial to improve the positioning accuracy [58]. In this article, we select 7 satellites in total to ensure positioning accuracy.

A. COMPARISON OF TWO SATELLITE SELECTING CRITERIA
In order to illustrate the effectiveness of the satellite selection criterion NWGDOP, we compared the satellite selection results based on the NWGDOP and GDOP criteria. There are 21 observable satellites including 10 GPS satellites and 11 Beidou satellites as Fig. 6 shows. And there are totally C 7 21 = 116288 subsets as 7 satellites are expected to be selected from 21 satellites. We randomly select 1500 satellite subsets from all subsets and compute their GDOP value and NWGDOP value. The positioning errors of the satellite subsets are computed and compared. Table 1 shows the positioning errors with GDOP and NWGDOP as the satellite selection criterion are obtained respectively. All satellite subsets are sorted according to GDOP or NWGDOP value. And we take the first 100 subsets with the minimum GDOP/NWGDOP value respectively. Fig. 7(a) shows the relationship between GDOP and RMS positioning error for 100 subsets with the minimum GDOP value. Similarly, Fig. 7(b) shows the relationship between NWGDOP and RMS positioning error for 100 subsets with the minimum NWGDOP value. As shown in Fig. 7, the RMS error doesn't change obviously as GDOP increases. And the RMS error decreases as NWGDOP decreases. The positioning error is more relevant to NWG-DOP in multi-system. The smaller the NWGDOP, the more accurate the positioning. So we can get the NWGDOP criterion to select optimal satellite subset in this article.

B. COMPARISON WITH THE RECENT STATE-OF-THE-ART METHODS
In order to verify the efficiency of this proposed method, the traditional satellite selecting method MGV, the basic ACO method, and two recent state-of-the-art methods including genetic algorithm (GA) and artificial neural network (ANN) are compared with the proposed algorithm ACO-PF algorithm in this section.
• MGV, in which the optimal satellite subset with the minimum GDOP value is selected from all possible subsets.
• ACO, in which the basic ACO algorithm is used to select the optimal satellite subset without any optimization. The criterion for satellite selection is NWGDOP.
• GA [10], in which GDOP measurement data is classified to GDOP threshold for classification. Thus a satellite subset with ideal GDOP can be selected. In this research we change the criterion GDOP to NWGDOP for a fair comparison.
• ANN [31], in which the GDOP classifier is employed for selecting one of the acceptable subsets. Similarly, for the sake of generality, GDOP is changed to NWGDOP in our compared experiment.
The comparison experiment compares the differences of the related methods under the condition of the same observation location and the same epoch.
First of all, to determine the optimal parameters of the comparison experiment, we made additional verified experiment on parameters setting of ACO-PF method. The five key parameters that have distinct effect on the proposed method are considered, which include the number of ants K , the heuristic factor α, the expectation heuristic factor β, the pheromone evaporation factor ρ and the pheromone intensity Q. The levels for five parameters are set as follows: four levels for K ∈ {10, 20, 50, 100}; four levels for α ∈ {0.5, 1, 2, 5} and β ∈ {0.5, 1, 2, 5}; four levels for ρ ∈ {0.3, 0.5, 0.7, 0.9} ; four levels for Q ∈ {1, 100, 1000, 10000} . The Taguchi's method adopts the orthogonal arrays so the orthogonal test table of L16(4 5 ) that contains only 16 experiments. Compared with the full-factorial analysis that needs 4 5 = 1024 experiments, Taguchi's method greatly reduced complexity. The specific experimental scheme is designed as Table 2 shows. We use the Taguchi's experimental design method in [59], which aims at finding the ACO-PF algorithm parameters with the minimum running time and the best robustness. The Taguchi's method determines the optimal combination of parameters by orthogonal design. It can comprehensively consider the relationship between parameters with less calculation, so as to find the reasonable combination of parameters [60]. As for ACO-PF method, the NWGDOP of satellite subset can indicate the effectiveness of the algorithm. And the number of iterations when the algorithm converges can be used as an indicator of calculation speed. Thus the two indicators including NWGDOP and convergence iteration number are set as the signal-to-noise ratio(S/N) in this experiment.
where x i , y i represent the output characteristics of the system, i.e., the minimum NWGDOP and convergence iteration number. And n is the number of levels divided in Taguchi's method, i.e., n =∈ {1, 2, 3, 4}.
The experimental results including minimum NWGDOP and the convergence iteration number of 16 sets of tests are also shown in Table 1. Put the minimum NWGDOP of each group of tests into (26) to get S/N NWGDOP , put the iteration number into (27) to get S/N iteration , and then add the two SNRs to get the required minimum SNR. As a result, we obtain the optimal combination of parameters : K = 10, α = 1, β = 5, ρ = 0.7, Q = 100.
As for the introduced two parameters, the iterations B of initial stage is generally no more than 60% of the number of convergence iterations. And the NWGDOP reference value r is usually a little lager than the minimum NWGDOP. According to the size of satellite selection problem in this article, we set all parameters in ACO-PF as Table 3 shows. After setting the value of the user-defined parameters, the comparison experiment between the proposed algorithm and other algorithms executes. Fig. 8 shows the changes of NWGDOP value in five compared methods during a fourhour period. Compared with other methods, the fluctuation of NWGDOP in ACO-PF method(the second subfigure) is steadier and almost all the values fall below 2.5. The NWG-DOP in the MGV method is obviously not the smallest, because the impact of noise on the satellite signal is not considered in MGV. Compared with the basic ACO algorithm, the change range of NWGDOP of ACO-PF is obviously narrowed, and the change trend tends to be stable. It can be seen that the improvements proposed in this article improve the performance of the algorithm. The NWGDOP performances of two recent state-of-the-art methods including GA and ANN are not superior to ACO-PF since they perform fluctuant. Fig. 9 shows the positioning errors by different methods and further reveals the characteristics of different methods. In Fig. 9, we take the mean value of RMS positioning error in an hour, and it lasts for 24 hours. It can be seen clearly that the RMS positioning error by ACO-PF method is smallest and the most stable among all compared methods. This further proves that the positioning accuracy of selected satellites by ACO-PF algorithm is higher.    Table 4 summarizes the experimental results including NWGDOP value, positioning error, and computational time of all related methods. As for average NWGDOP and RMS positioning error, the ACO-PF performs better than other methods. The satellite subset selected by ACO-PF algorithm has high accuracy. The CPU time listed in Table 4 is also the average runtime consumption of five methods in one cycle over 24-hour experiment. We can see that the GAbased algorithm in the experiment of selecting 7 satellites from 21 satellites has the lowest time consumption. Due to the pattern training time, the ANN algorithm need more time for execute an iteration. Compared to these two methods, the time consumption of the ACO-PF algorithm has not increased significantly and it has been significantly reduced compared to the basic ACO algorithm. Therefore, the ACO-PF satellite selection algorithm proposed in this article has high positioning accuracy, stability and acceptable runtime consumption.

C. ELABORATION ON THE EFFECT OF TWO IMPROVEMENTS
Two variants of ACO-PF including ACO-PF-1 and ACO-PF-2 are constructed by incorporating only one of the two improvements to clearly show the impact of each improvement.
• ACO-PF-1, which is the variant of ACO-PF by incorporating a polarized feedback. Through the comparison between ACO and ACO-PF-1, we can get the effect of polarized feedback.
• ACO-PF-2, which is the variant of ACO-PF by incorporating a perturbation operator. Relatively, ACO and ACO-PF-2 are compared to get the effect of a perturbation operator on the algorithm.
The comparison experiment between ACO, ACO-PF-1, ACO-PF-2 and ACO-PF mainly explores the effect of each specific improvement on algorithm's convergence performance and searching result. Fig. 10 and Fig. 11 show the convergence performance comparison of related methods over a large number of satellite selection experiments. Each method has been carried out in 200 times experiments. In Fig. 10 and Fig. 11, each point represents an independent experimental convergence result, and the corresponding abscissa value indicates the iteration number when the program converges; the ordinate value indicates the NWGDOP value of the solution, that is, the NWGDOP value of the satellite combination.
Comparing the distribution of the searching results of ACO-PF-1 method (pink triangle points) and ACO method (dark blue square points) in Fig. 10, we can find that the results of ACO-PF-1 are distributed on the left side of the area. According to the abscissa value, the numbers of convergence iteration are mostly within 80, and the average convergence iteration is lower than that of the ACO method. The above analysis shows that the polarized feedback added in the basic ACO method effectively improves the convergence speed of the algorithm. Similarly, the results of ACO-PF-2 (light blue dots) are distributed on the lower side of the area. The NWGDOP values of ACO-PF-2 method are almost in the range between 2.20 and 2.24, which are far smaller and more concentrated than the results of the ACO method (2.23 2.31). It shows that the another added improvement, i.e., a disturbance factor, improved the searching ability of the basic ACO algorithm to search for better results.
In Fig. 11, we compared the ACO algorithm and the ACO-PF algorithm in the same graph. It can be seen that the distribution of 200 points representing ACO algorithm is VOLUME 8, 2020  dispersed, and 200 points representing ACO-PF algorithm are mainly distributed in the lower left area of the coordinate system. This shows that the ACO-PF method can achieve algorithm convergence within a small number of iterations (within 80 generations), and the performance of the solution is better (the NWGDOP value is smaller). It can be seen that the polarized feedback strategy and the designed perturbation operator proposed in this article effectively improve the convergence speed and the global searching ability of the algorithm, making ACO-PF an efficient satellite selection method.
Additionally, we made statistics about the relationship of NWGDOP value and iteration index over a large number of experimental results. We conducted a total of 50 sets of experiments, each set of experiments were conducted 200 random realizations similar to that in Fig. 10 and Fig. 11. We record the average NWGDOP value of 50 results, forming Fig. 12. The maximum and minimum value of NWGDOP at some iterations of the 200 times experiments are also marked. As illustrated in Fig. 12, it is clearly that ACO-PF-1 method (pink curve) converges faster than ACO method (dark blue curve) and the optimal NWGDOP reduces to a certain extent even if it is not the smallest. And the optimal NWGDOP value of ACO-PF-2 method (light blue curve) is significantly reduced compared with ACO method, which is nearly to the NWGDOP of MGV method. The performance of the two variants of ACO indicates the corresponding improvement effect, that is, a polarized feedback effectively improves the convergence speed of the algorithm and a disturbance factor improves the searching ability of the basic ACO algorithm.
Through the separate comparison of ACO method and ACO-PF method, we find that ACO-PF method converges speedily within 50 generations than the original ACO method which converges around 80 generations. We can obtain the convergence speed of ACO-PF is more 50% faster than the original method. When the two methods converge completely, the NWGDOP value of ACO-PF method is 0.065 smaller than the basic ACO shown in Fig. 12. The inset figure of Fig. 12 shows the distribution of average NWGDOP value for 200 experimental realizations when iteration is 120 (in dotted ellipse). Obviously, NWGDOP values distribution of the proposed method in this article is more concentrated and most of them distribute within the average optimal range. Fig. 12 shows the fast convergence speed and the stable searching performance of ACO-PF method relative to ACO method.

VI. CONCLUSION
In this article, we discussed the navigation satellites selection problem, which is presented as a combinatorial optimization problem. A navigation satellites selecting method based on ant colony algorithm with polarized feedback has been proposed to select the optimal or near optimal satellites subset. And the NWGDOP is defined as the criterion of satellite selection. The NWGDOP considers not only the geometry of the satellite combination but also the satellite's elevation and SNR which influence positioning accuracy. The perturbation operator designed in the selection probability in algorithm process can improve the global searching ability of the algorithm. At the same time, the proposed polarized feedback mechanism can greatly improve the convergence speed of the algorithm. Through extensive experiments, we demonstrate that the proposed algorithm has high convergence and powerful ability of global searching compared to the original ACO algorithm, GA as well as ANN algorithm. In terms of convergence, ACO-PF outperforms the original method by up to 50%. When the algorithm converge completely, the NWG-DOP value of ACO-PF is always 0.065 smaller than the original method. That demonstrates that the global searching ability of ACO-PF is much closer to the optimal solution. The proposed improved algorithm is applied to the selection of navigation satellites in this article and outperforms the related methods in navigation satellites selection so it will have broad application prospects. The theoretical provement about the convergence performance of the proposed method is also expected and left for further investigation. He is currently an Associate Professor with the School of Computer and Information, Hefei University of Technology. His current research interests include multi-agent systems, evolutionary computation, and search-based software engineering.
YONGTANG YU received the Ph.D. degree from the Xi'an University of Architecture and Technology, Xi'an, China, in 2020.
Since 2016, he has been an Associate Professor with China Jikan Research Institute of Engineering Investigations and Design Company Ltd., Xi'an. His research interests include geosensor networks and geotechnical data processing. He is an individual member of the CCES.
JIWEN ZHANG received the M.Eng. degree from the Xi'an University of Technology, Xi'an, China, in 2004.
Since 2010, he has been a Professor with China Jikan Research Institute of Engineering Investigations and Design Company Ltd., Xi'an. His research interests include sensing technologies and information systems in geotechnical engineering. He is an individual member of the CCES and CECS. VOLUME 8, 2020