Robust Estimators in Geodetic Networks Based on a New Metaheuristic: Independent Vortices Search

Geodetic networks provide accurate three-dimensional control points for mapping activities, geoinformation, and infrastructure works. Accurate computation and adjustment are necessary, as all data collection is vulnerable to outliers. Applying a Least Squares (LS) process can lead to inaccuracy over many points in such conditions. Robust Estimator (RE) methods are less sensitive to outliers and provide an alternative to conventional LS. To solve the RE functions, we propose a new metaheuristic (MH), based on the Vortex Search (IVS) algorithm, along with a novel search space definition scheme. Numerous scenarios for a Global Navigation Satellite Systems (GNSS)-based network are generated to compare and analyze the behavior of several known REs. A classic iterative RE and an LS process are also tested for comparison. We analyze the median and trim position of several estimators, in order to verify their impact on the estimates. The tests show that IVS performs better than the original algorithm; therefore, we adopted it in all subsequent RE computations. Regarding network adjustments, outcomes in the parameter estimation show that REs achieve better results in large-scale outliers’ scenarios. For detection, both LS and REs identify most outliers in schemes with large outliers.


Introduction
In 2000, the U.S. government turned off the selective availability of the Global Positioning System (GPS), making it more responsive to civil and commercial use. Since then, the industry of GNSS-based technologies (Global Navigation Satellite Systems, which includes other satellite constellations) has grown significantly. Surveyors, researchers, and most civilians use products developed from this technology. At present, GNSS receivers can be found in many devices, from precise measuring instruments to cellphones and cars. In the field of geodesy, GNSS has facilitated the implementation and quality control of high-precision and accurate networks. Before the advent of GNSS, establishing a geodetic network required a direct intervisibility between the points, as in triangulation and traverse surveys [1]. Now, point co-ordinates can be defined by GNSS signals coming from orbital space, which only requires a good coverage of satellites during data collection. This has brought great flexibility in network design, and the achievable accuracy has also helped the densification of classical first-order geodetic networks [2].
Most geodetic networks have been established using GNSS signals, which are available all over the Earth's surface. They are defined by several materialized points with high precision and their co-ordinates serve as a basis in a wide variety of applications. Some examples are: support and basic control for surveying and mapping projects; implementation and maintenance of infrastructure works; monitoring of structural deformations; land cadastre and management; and monitoring of geodynamic events, such as earthquakes, landslides, and volcanoes. The study of Mt. Etna is a good example; two GNSS-based networks found a motion pattern of the active Nizzeti faults [3]. In addition, a geodetic network design was applied to find optimal ground locations for interferometric synthetic aperture radar (InSAR) devices [4].
To establish a geodetic network using GNSS, the receivers need to be placed at benchmarks and record satellite tracking data over a given period. Then, baseline vectors are determined between pairs of receivers, based on ranging and phase observations to the satellites. Once the vectors are computed, they become the new observations of the network and need to be adjusted to get the final point co-ordinates [1].
Although most of the GNSS observation process does not require much direct human influence, the vectors may not be free of blunders. There are sources of error that may cause significant deviations (e.g., multipath signal propagation error, cycle slips, and ionospheric anomalies). The electromagnetic wave assumed to propagate along the line of sight between the satellite and the receiver might be reflected or scattered by obstructions before reaching the instrument [5]. Other error sources are more mundane, such as mistakes in measuring the height of the antenna above the marker. Although random errors are inherent to observations, outliers should be detected, identified, and adapted for better determination of the co-ordinates [6]. An outlier is better not to be used (or not used as it is) in an adjustment process because it has a high probability of being caused by a gross error [7]. In geodesy, outliers are mostly produced by gross errors, and gross errors most often lead to outliers.
As outliers in observations affect the accuracy, those errors need to be identified or minimized [8]. In standard geodetic adjustments, a least squares (LS) process is often applied, as it is the best linear unbiased estimator, assuming that no outliers and/or systematic errors are present in observations [1,9,10]. If the data are contaminated, however, such an estimation will lead to biased parameters [11].
Alternatively, the adjustment can be based on Robust Estimators (REs), which are less sensitive to outliers. If some observations contain blunders or even systematic errors, the REs will be insensitive to those non-random errors when estimating the parameters [11]. REs have a wide range of applications; for example, in [12], a robust parameter-estimation method for a mixture model working with the weights of samples is presented. Furthermore, applications of non-Gaussian distributions on multipass SAR Interferometry have been presented in [13], and new methods for geodetic observations have been presented in [14]. Therefore, they need to be investigated for a better understanding of their capabilities and limitations in geodetic networks.
The minimizing objective functions associated with REs are not linear and, therefore, iterative processes or smart techniques are required to solve them. One approach is to combine metaheuristic algorithms (MHs), to optimize the objective function. This strategy has been applied in several studies in geodesy, as it may lead to better results than classical methods [15][16][17][18][19][20][21]. In MH research, the particle swarm optimization (PSO) [22] has been widely applied, followed by the artificial bee colony [23][24][25] and ant colony optimization methods [26], more recently [27]. Another solution is to apply a data-snooping procedure [28,29], a statistical test which takes place after the LS computation. However, studying statistical tests is not the purpose of this paper. Researchers have also applied MHs in others areas of geoscience, such as remote sensing [30,31].
For the adjustment of geodetic networks, MHs compute the unknown parameters (point co-ordinates) within a pre-defined search space and check them by evaluating the objective function for the chosen RE. The estimate is then used for the goodness-of-fitting evaluation. The process is repeated, following the exploration strategy adopted by the MH, until either an acceptable solution or the computational limit is reached [20].
Applying MHs in geodetic networks adjustment has not been widely explored. Most studies have been limited to working with only one or two REs, often in simple cases. For example, they may not have considered the generation of multiple error scenarios [18,19,21], or may have omitted the consideration of random errors [20]. It is important to explore adjustment solutions where random errors can be simulated, getting as close to a real situation as possible. Still, several error scenarios need to be tested with a variety of REs to obtain a representative set and measure the reliability of the strategy.
Furthermore, an MH needs to present stable results when minimizing the RE functions and present good strategies to avoid local minima.
Therefore, this work has two main goals: First, developing better MHs by presenting the Independent Vortices Search (IVS) based on the Vortex Search method [32], which brings advances to mitigate the limitations found in the original algorithm [33]. We compare its performance with the original method, other modifications, and a competing algorithm.
The second goal is improving quality control in geodetic networks by combining smart and robust methods. By separating the results into three scenarios (with no outliers, small outliers, and large outliers), we classify the REs according to the demand and situation.
We generated observations from the official co-ordinates of a station and standard deviations from GNSS signal processing. Several error scenarios were simulated in the observations, which are tested by all REs, and the results are compared with the conventional LS. In addition, different median and trim positions in residual vectors of some REs are investigated. We analyzed the results by: (1) the identification outliers, analyzing the residuals vector, and (2) by the estimated parameters, comparing them with the official co-ordinates of the network.
Additionally, for problems based on the Gauss-Markov model, we present a new proposal for the search space definition when using MHs. Furthermore, a comprehensive analysis of the results achieved with different trim and median position in some REs is given, showing which equation is best in each situation. This paper is organized into six sections: In Section 2, we present a theoretical overview of the LS and RE methods applied in the experiments. Section 3 presents the first contribution of the study, showing the Independent Vortices Search algorithm. In Section 4, we present details of the experimental setup, how we analyzed the quality of the solutions, and the implementation of the REs and IVS. Section 5 presents the second contribution of this work, discussing the results obtained with adjustments to the GNSS network. We compare the results of the trim and median position of some REs and the outcomes, as well as discussing the various error scenarios. Section 6 brings the final considerations about the research and suggests issues to explore in future research.

Classical and Robust Approaches
The adjustment of geodetic networks is generally conducted by employing a Least Squares (LS) method, which is the best linear unbiased estimator when only random errors are present in the observations [9]. When this assumption does not comply, the estimation fails, distributing errors over many parameters. Other techniques can be applied when one suspects the presence of outliers, such as robust estimators (REs). We present a brief explanation of the methods covered in this study below.

Least Squares Method
The adjustment of geodetic networks applying the LS method is given by the Gauss-Markov Model [1]: where A is the design matrix of the partial derivatives, x is the parameter correction vector (for point co-ordinates) of the linearized model, l is the observations vector minus the calculated data, and v is the residuals vector. The calculation of A results in a matrix formed by the values 0, 1, and −1, based on the partial derivatives of the observation equations: representing the distance between the approximate parameters on each axis. The solution of the conventional LS adjustment is given by [9]: where T represents the matrix transpose and −1 represents matrix inversion. The quantity W is the weight matrix of the observations [1]: where Σ is the covariance matrix of the observations and σ 2 is the a priori variance factor. The final solution is given by Equation (5): where x is the parameters' correction for the initial approximate parameter vector (x 0 ), andx the vector of adjusted parameters containing the final co-ordinates of the network points.

Least Trimmed Squares
The Least Trimmed Squares (LTS) robust estimator consists of summing only the h smallest squares of the residuals (Equation (1)), suppressing the larger ones [34]. The estimator is given by: According to ( [38] (p. 132)), h is defined by: where n equals the number of observations and p is the number of parameters, which is the number of co-ordinates to be estimated.

Least Median of Squares
The Least Median of Squares (LMS) replaces the squared residuals summation of LS with the median square residual [34], given by: where the median position (med h ) is defined by Equation (7).

Least Trimmed Squares with Redundancy Constraint (LTS-RC)
The LTS-RC extends LTS by adding a constraint to the summation of residuals. This estimator ensures that each parameter remains associated with at least two observations represented in the summation by their corresponding residuals. This measure guarantees observation redundancy [20].
For the adjustment of geodetic networks established only with GNSS baseline vectors as observations, the algorithm maintains at least two residuals associated with each parameter, guaranteeing redundancy. In estimation problems where a parameter depends on more than one type of observation, however, the algorithm must be adapted [39]. This is to ensure that each different type of observation is present at least twice.

Sign-Constrained Robust Least Squares
Unlike other REs presented so far, the Sign-Constrained Robust Least Squares (SRLS) has multiple optimization functions [35]: where the diagonal elements of W R are defined as with C being a positive constant defined in the range [1][2], according to the estimated number of outliers in the observations. For cases with 1% of outliers, C = 2.0; for 10%, C = 1.0. The quantityσ is a robust estimate of the standard deviation of the unit weight, in practice given by [40]: where med √ w i |v i | is the median absolute values of studentized residuals. The number r W R is the rank of W R . However, due to the presence of covariance values in the weight matrix of GNSS networks, the matrix W R can be given by [18,37]: where diag w R 1 , . . . , w R n are the elements calculated by Equation (10), distributed along a diagonal matrix.
In addition to the minimization of Equation (9) and the constraint for W R , the SRLS is also subject to the following summation [35]:

Least Trimmed Absolute Deviations
The Least Trimmed Absolute deviations (LTA) have a formulation similar to LTS. However, the squares of the residuals are replaced by their absolute values [36]. A summation of the h smallest absolute residuals is carried out, ignoring the largest ones: where h is defined as in LTS, using Equation (7).

Iteratively Reweighted Least Squares
The Iteratively Reweighted Least Squares (IRLS) method represents a class of REs that works on the reweighting of W to estimate the parameters. This type of estimator has its own iterative methods to obtain a robust solution, thus not needing any external numerical method (such as MH).
The iterative method chosen was proposed in [37] and has presented good results [20,37,41].
In [37], the author presented two weight adjustment strategies. Strategy I held better results and, therefore, is adopted in this study: where τ i is given by with c i being a unit vector, only filled with a 1 in the ith position, and withσ given by Equation (11). The quantity R is the redundancy matrix [42] where I is an identity matrix. Strategy I is, then, computed as follows [37]: 1. Initialize the counter i = 0 and calculate the initial solution of the parametersx (0) using Equation (3). 2. Increment i = i + 1 and implement the weight adjustment strategy I using Equation (15). 3. Construct the equivalent weight matrix W R using Equation (12) by applying the calculated values to w i . 4. Stop the iteration if the difference x (i) −x (i−1) is less than a threshold value, or if a maximum number of iterations has been reached. Otherwise, return to step 2 and start a new iteration.

LMS Median Position and LTS/LTS-RC Summation Limit
As seen above, the LMS has a median position defined by h by applying Equation (7). Likewise, the trim limits (residuals for summation) for LTS and LTS-RC are defined, based on the same equation.
Occasionally, a variation of this equation is used [19,21,43], replacing Equation (7) by Equation (18) [38], below: where n is the number of observations. It is clear that the only difference, regarding Equation (7), is the absence of the variable p (the number of parameters). By defining h as above, more residuals are trimmed in LTS and LTS-RC, and the true median position of the residuals vector becomes the median value for LMS.
We can also think about other possible values for h. In the case of LTS and LTS-RC, a lower limit considers fewer residuals in the estimation (i.e., it will ignore more values in the residuals vector). For the LMS, the limit defines a single residual position to be taken into consideration. With a greater limit, the residual location considered the center will be further from the median, toward the larger values.
Thus, besides testing and analyzing the REs presented in this study, we also will verify the behavior of both Equations (7) and (18) for LMS, LTS, and LTS-RC estimators.

Independent Vortices Search
All of the REs (except for IRLS) require a numerical optimization method. The MH has to explore values for the parameters which generate the residuals from the Gauss-Markov model. The objective function of the RE is then evaluated and used for the computation of the solution fitness. Then, the MH generates new candidate solutions, continuing the process until it reaches a pre-defined limit [20].
For the numerical optimization, we developed the Independent Vortices Search (IVS), based on the Vortex Search algorithm (VS) [32] and its modification, the Modified Vortex Search algorithm (MVS) [33]. While VS works with one vortex, MVS can explore many vortices at once, exchanging information. We will present more details below.

Vortex Search Algorithm
The Vortex Search algorithm (VS) is a single solution MH for numerical function optimization. It generates new candidates around the current best result, moving along the search space when achieving a better solution. The radius defines the limit for generating new candidate solutions and decreases over cycles [32].
VS and MVS define the initial radius (radii) of the center (centers) by using Equation (19), covering the whole search space. The algorithms obtain their first solution by applying Equation (20), causing the search to begin at the center of the space: where r 0j is the initial radius of the jth parameter, x min,j and x max,j are the respective minimum and maximum limits of jth parameter, and µ 0j is the center of the vortex and the initial solution for the parameter j.
In each cycle, the radius decreases to limit the generation space for new candidate solutions. In the later cycles, the VS produces a fine adjustment as the current candidates get closer to the best solution. The radius decrease is obtained by Equation (21), and provides satisfactory control of the exploration and investigation [32]: where r c is the radius size in cycle c, R is a constant (set at 0.1) that controls the resolution of the algorithm's search, and γ (R, a c ) is the incomplete inverse gamma function, given by Equation (22): where a c is a parameter that defines the shape for each cycle, as defined by Equation (23): The equation of the shape parameter uses the cycle count c and the maximum number of iterations C max .
The vector of candidate solutions s is randomly generated at the center µ using the Multivariate Normal Distribution (MND) with the standard deviation vector of r. Figure 1 shows an example of the radius convergence by applying Equation (21) and the candidate generation limits by the MND. The algorithm checks that the candidate values are within the search space by applying Equation (24), where is a random real value between [0, 1]: For each candidate solution, a fitness function tests the result: If a candidate produces a better outcome than the current one, VS moves the center (µ) to the new solution and discards the remaining candidates. Otherwise, the algorithm keeps the best obtained result and produces new solutions.
Both VS and MVS use all the above equations. However, by using several vortices, MVS adjusts the positions of its centers to each cycle, based on the vortex, with a better result. This separates positioning the center µ l for the generation of new solutions, and the best solution s l of a vortex l. For this, MVS applies Equation (25): where µ l is the new position of the center l, s l is the best solution of the center l, and s best is the result of the vortex with the best solution. At each cycle, all centers (except for s best ) gain new positions using the above equation [32]. As has been shown, VS has no crossover nor mutations between candidate solutions, unlike Genetic Algorithms (GAs) and the Artificial Bee Colony (ABC) method. At each cycle, it saves the best solution-either the current best or a new one-for the next cycle, discarding the remaining ones. This way, VS has no need to use individual selection strategies, a mutation rate, or other strategies present in other MHs. This helps to keep the configuration simple, reducing human interference in searching for a solution to a problem.
Another interesting characteristic of VS is the transition between global and local search. In many MHs, exploration and investigation occur from the beginning to the end of the execution. This makes it possible to find a good enough solution, even in the first few cycles. This way, they can interrupt the execution if so desired. Otherwise, the global search continues through to the last iterations, attempting to reduce the risk of local minima. In VS, the radius size is reduced over the cycles using an inverse incomplete gamma function (Equation (22)) and considering the ratio of cycles performed (Equation (21)). Thus, exploration and investigation occur gradually, and so are not present throughout the entire execution, yet still providing an adequate balance. The advantage of this strategy is the fine-tuning of the best solution found in the last cycles.
However, by analyzing the VS strategy, we can conclude that the transition between exploration and investigation has two drawbacks. First, it has the local minimum risk: There is no possibility of achieving the global solution in an execution if it limits the radius to a local minimum area. The radius of the vortex becomes small and the candidates get stuck in a local search limitation. The second inconvenience is identifying the dispensability of the processing during the execution. When defining a limit of execution for example, by the number of cycles, it executes them in its totality because only in the later cycles does VS carry out the investigation stage.
MVS tries to minimize the local problem by adding several vortices to VS. However, applying Equation (25) does not mean the vortex with the so-far best solution will bring the other centers closer to a global minimum.
With these characteristics in mind, we based the proposed IVS method on maintaining solution fine-tuning while trying to circumvent the local minimum problem, all without increasing the computational cost.

Characteristics of the IVS
The IVS starts from the MVS, aiming to make it simpler and more efficient. We propose three strategies through which IVS differs from MVS.
The first major difference is the complete elimination of Equation (25), which changes the vortex center positions, based on the vortex with the best result. By doing so, each vortex can proceed toward a minimum, whether local or global. This is intended to minimize the local solutions, as we can add more vortices to explore solutions independently. Figure 2 presents a hypothetical situation in which the vortices follow their paths towards the minima.  The second feature of IVS is the replacement of Equation (20) for each additional vortex. The first vortex has an initial solution defined by applying Equation (20), whereas the remaining, by Equation (26), are randomly distributed in the search space: where µ l 0j is the initial solution for the jth parameter of the vortex l, with l > 0, x min,j , and x max,j are the respective minimum and maximum limits of the jth parameter , and is a random real value in the interval [0, 1].
The third modification considers maintaining the amount of candidates per vortex when L > 1 (where L is the total number of vortices). In the MVS proposal [33], the experiments presented five centers and divided the number of candidates between them. IVS, however, fixes the quantity of candidates for each vortex. To compensate for the computational increase, IVS decreases the number of cycles performed to keep the total Fitness Evaluations (FEs) the same.
It is still possible to verify that the first two strategies have the possibility of a higher diversity of the candidates generated. By removing Equation (25), each vortex is allowed to go ahead with its search, not dragging them close to a vortex, with better overall results at the moment. The second strategy (Equation (26)) allows for greater diversification while generating the first results. This is because the MND, which uses the center of the vortex as the mean for candidate generation, generates more solutions near the center. With the centers of vortices more widely distributed, IVS should not concentrate candidates in only a part of the search space.

A New Search Space Definition
For the definition of the search space in the adjustment of geodetic networks, we defined the limits of each parameter j by x min,j and x max,j . By using the parameter vector x obtained from conventional LS computation, we set the exploration area according to Equation (27): where max i |v i | is the greatest absolute residual from the LS solution, max i diag (R e ) i is the greatest absorption fraction of an error and 1.15 indicates a 15% safety margin. Additionally, R e is defined by Equation (28) [44]: where Σv is the covariance matrix of adjusted residuals from the LS solution. Equation (27) presents a novel method to define the search space. It can be applied to any problem susceptible to outliers that use the Gauss-Markov model to generate residuals and MH to explore the solution. It is not limited to geodetic networks and can be used along with REs or other techniques-for example, in hyperspectral image data [45,46] or 3D Point Clouds [47]. Applying this strategy to 900 scenarios, it did not omit any solution in the IVS exploration limit.

The GNSS Network and Simulations
To perform the tests, we built a geodetic network from six stations in the Brazilian Network of Continuous Monitoring of Global Navigation Satellite Systems (RBMC) (Figure 3). The observation vectors were formed from the official co-ordinates of the points. We processed 6 h of data from the stations to obtain the covariance matrix (Σ y ). This allowed us to compute the weight matrix (W) and generate random errors for the observations. The vector l was set up by adding the official values and the random errors generated from the covariance matrix. The network was structured without repeated observations and with independent baselines. This gives six points (one control point) and 13 vectors. As they are three-dimensional, there were a total of 15 parameters and 39 observations.
To test the RE performances, we built nine packages with 100 error scenarios each, as presented in Table 1. Package 01 had no outlier, only random errors in all observations ranging from [0σ, 3σ) using a normal distribution. Trimming random errors at 3σ should not have affected the results. Following the normal distribution, the occurrence of values greater than 3σ was only 0.27%. All estimators were tested in this condition.
The other packages contained at least one outlier per scenario. Outliers were also calculated from the standard deviation of each observation. The packages separated small error scenarios (with magnitude between [3σ, 6σ)) from large error scenarios ([6σ, 12σ]). Outliers were randomly distributed over l. Packages with four outliers represented scenarios with approximately 10% contaminated observations. In the most critical cases, the network should be resistant to two simultaneous outliers between the same vertices, as we had at least four vectors connecting each station.

Analysis and Validation
In this work, we analyzed the quality of the solutions of each RE and the conventional LS method in three ways. First, we classified the solution regarding whether the outliers were detected or not. Second, we quantified detected, false positive, and unidentified outliers. Third, the numerical solution of the estimator was compared to the true and known co-ordinates of the points. The two first analyses were obtained using the residuals vector v, while the numerical comparison was extracted from the estimated co-ordinates vectorx.
We divide the classification into six classes, represented by capital letters from "A" to "F": We considered false positive to be a value bigger than 3σ in the residuals vector in a position where no outlier was inserted, as we truncated the random errors by up to three sigma.
Comparison of the estimated solution with the true co-ordinates was carried out in order verify the impact of the errors in the final solution of the adjusted parameters (x). This allows us to consider, for example, whether an estimator was better for outlier detection or parameter estimation.

IVS Parameters
The parameters of IVS were the same for all scenarios and REs. A total of 50 candidates per vortex and 40 vortices were used, as this configuration gave the best results in our experiments. To limit the algorithm execution time, we defined the amount of fitness evaluations (FE) to 25,000,000.
The fitness evaluation of the RE solutions was conducted by Equation (29) [20]: where f is the result of the objective function.

Robust Estimators Adaptation
In LS, the weight matrix of the observations (W) plays an important role for the estimation; therefore, it is also necessary to consider it here. In GNSS networks specifically, observations are correlated, rendering the problem even more complex [41]. These observations form a weight matrix with some negative values of covariance. Hence, in LTS, the summation of the residuals has to be performed after extracting the absolute values of the weighted squared residuals, in order to eliminate the influence of the signs in the ranking. This results in the following Equation (30): where w i is defined by with σ v representing a permutation of the indices i = 1, . . . , n, such that 0 ≤ w 1 ≤ w 2 ≤ · · · ≤ w n . The quantity v 2 represents the vector that has, as components, the squares of the residuals in v.
The inner product is, then, given by where W i = (W i,1 , W i,2 , . . . , W i,n ) denote the respective rows of W. Following this idea, for LMS, the mathematical model already adapted to GNSS networks is given by: where w i is given by Equation (31), and the median position (med h ) is defined by Equation (7). Likewise, the LTA adapted for GNSS networks is given by: with z i following the same idea as the LTS, except that v is not squared: For the SRLS, Equation (13) presents a constraint that needs to be implemented. In this study, we adopted a penalty function that multiplies the result of Equation (9) given by Equation (36): where f p is the penalized result for the SRLS estimation. If Equation (13) equals 0, f is multiplied by 1, and the result remains unchanged. Otherwise, it will be penalized. In the sign function, absolute values smaller than 1 × 10 −9 were considered as 0.

Flowchart
The development and implementation of the adjustment routine using REs and IVS followed the flowchart shown in Figure 4.

Results and Discussion
For better comprehension, this section is divided into three subsections: Section 5.1 presents the optimization with IVS and comparison with the performances of other MHs. Section 5.2 shows the estimated impact of the trim limit and median position in some REs, when applied to the geodetic network adjustment. Section 5.3 presents the results for the network adjustments, applying IVS and testing all the REs (including the trim limit and median position variations) in 900 error scenarios. The adjustments were organized in three topics: no outlier case, small magnitude outliers, and large magnitude outliers.

Performance and Discussion of IVS
We conducted several tests to compare the solutions obtained with IVS to solutions obtained by other MHs. For a better analysis, we generated four scenarios of errors in a GNSS network. Each scenario contained three simultaneous outliers, ranging from [3σ, 6σ), in 39 observations (about 7.8% contamination). Using the LMS estimator, the MHs executed the scenarios twenty-five times. This allowed for analyzing the stability and quality of the results. More details of the network can be found in the next section.
Besides the MHs already presented, we also tested the Hybrid Vortex Search algorithm (HVS). The HVS works with the combination of ABC and VS, trying to take advantage of the best characteristics of each strategy. The reader can find more details in [49]. Thus, the algorithms tested were: • Artificial Bee Colony (ABC) [23,24]; • Vortex Search algorithm (VS) [32]; • Modified Vortex Search algorithm (MVS) [33]; • Hybrid Vortex Search algorithm (HVS) [49]; and, • Independent Vortices Search (IVS).
The configurations of the MHs obeyed the following parameters: 5,000,000 FE as processing limit, 50 candidates per center for vortex algorithms, and 50 bees in ABC (25 food sources). In addition, we tested the MVS and IVS with five vortices, as presented in [33], and 40 vortices. Preliminary tests showed better results from 20 to 60 vortices, so we adopted an intermediate value of 40. It is important to note that cycle counting was not used as a configuration parameter or as a limit for the MH executions. For more objective comparisons between MHs, we should not use the cycle count as a stopping criterion. As each MH has different strategies in its execution, it can use a different quantity of FE in each cycle. Using a fixed number of cycles as a parameter to limit the execution can cause a large variation in the total amount of FE for each MH. This would favor MHs that make more FE per cycle, since they have more opportunities to test their solutions [50]. We adapted the algorithms to use the number of FE as the limit for suspension. This allowed for a more adequate performance comparison. To better situate the reader, executing 5,000,000 FE was equal to 100,000 cycles in ABC with 50 bees (25 food sources); 100,000 cycles in VS with 50 candidates; and 20,000 cycles in IVS with five vortices and 50 candidates per center.
The MVS experiments performed in [33] divided the quantity of candidates by the number of vortices. Thus, in executions with 50 candidates and five vortices, MVS generated only 10 candidates per vortex every cycle. As this work uses the number of FE as the execution limit of the MHs, each vortex generated a number of candidates established by the number of candidates parameter. For a run with 50 candidates and five vortices, the algorithm generated 50 candidates per vortex each cycle. This characteristic is standard in IVS.
For better representation of results, we present them in box plots. Each point represents the result of one run. The number after "c" in brackets (e.g., "c [5]" and "c [40]") is the amount of vortices, where applicable. Figures 5-8 show the results for Scenarios 01, 02, 03, and 04, respectively.    By analyzing the charts, we see that IVS performed better than other MHs. The configuration with 40 vortices overcame the alternative with five vortices, although both options showed good results. IVS achieved the lowest minimum and presented little variation. The HVS strategy did not improve much on the solutions, as they were similiar to those found by the original VS. ABC showed lower variation in the results when compared to the already known MHs. This shows a certain stability in the solutions, surpassed only by IVS.
It is also possible to note the horizontal alignments of points (solutions) in the charts. This was common to several MHs and identifies the local minimums in which the MHs got stuck in, on some runs. In Figure 7, it is possible to notice an alignment between the values 2, 4 and 2, 6 for MVS and HVS. The same diagram shows another alignment between 1.8 and 2.0 for all MHs, which coincided with the bottom of the boxes for the VS, MVS, and HVS solutions.
The MH proposed in this work overcame all other tested MH in both the mean and variation of the solution. It needed no extra computation, producing the first scientific contribution of this work. We applied the IVS configuration with 40 vortices in the following experiments to optimize the RE functions.

Comparing LTS/LTS-RC Trim Limit and LMS Median Position
We performed tests with the two values for the limits of LTS and LTS-RC and for the median position of LMS. The lower limit, h = (n + 1) /2, will be shown as h n2 and the upper limit, h = (n + p + 1) /2, as h np2 . Figure 9 shows that the adjustments with the limit of h n2 had worse results in the classification for LMS and LTS. Classifications of type "A" were reduced by adopting the expected threshold, as compared to more robust, increased "B" classifications. For these estimators, the lower limit maintained outliers identification, but pointed out more false positives. The results for the LTS-RC present a slight improvement to the limit of h n2 . Comparing the results of LTS and LTS-RC, we notice that the results presented almost the same classifications, whereas the variations with h n2 presented a significant difference. This is because the larger limit (h np2 ) removed fewer residuals, hardly breaking the redundancy. Furthermore, since the network had several observations among its three-dimensional points (although not repeated), it was hard for any parameter to remain without at least two residuals (not breaking redundancy). Using the limit of h n2 , the classic LTS more easily broke the redundancy in the residuals trim, whereas the LTS-RC achieved more satisfactory results, retaining the redundancy of the network. Figure 10 confirms a gain of false positives with the limit of h n2 for LMS and LTS, and stability with a slight improvement for LTS-RC. There were no significant differences in outlier identification for the different limits, while, for undetected errors, only the LTS showed a 6.6% worsening.  We also compared the distance of the co-ordinates from their true values. Figure 11 displays the mean distance difference of the co-ordinates in each scenario package, between the limits h n2 and h np2 , for each RE. The mean differences were negligible in most situations, being below 1 mm with either h n2 or h np2 . LTS produced inferior results with h n2 , whereas the LTS-RC reduced the variation in two cases, with one and four simultaneous outliers. The redundancy constraint has proven to be interesting for any problems where the number of observations per parameter is more limited and/or the trim limit is smaller. To the best of our knowledge, this is the first time the h value has been analyzed for a geodetic application. As the measurements of these networks do not contain data related to all the parameters, the redundancy constraint also played an important role. Besides this specific treatment, the h value can be considered in any studies that wish to apply LMS, LTS, or LTS-RC.

Results in Network Adjustments
For better analysis of the behavior of each RE, we divided the results into three topics: the no outlier case (Package 01); scenarios with small outliers (Packages 02, 04, 06, and 08); and scenarios with outliers of great magnitude (Packages 03, 05, 07, and 09).  This means that all REs presented false positives in some scenarios. Figure 13  To check the solutions of each estimator, we elaborated the box-plot chart shown in Figure 14. The LS exceeded all other estimators, showing a lower variance and a lower variation between the scenarios. As expected, the LS estimation was superior to any other RE, when tested in the conditions of no outliers. In addition, LS required no MH to estimate the parameters and presented low computational cost, when compared with estimators that need MHs. It is important to point out that the points will not reach zero, due to the random errors present in all observation vectors. Both for outlier identification and co-ordinate estimation, LS proved to be the best method for scenarios with no outliers. The high number of false positives in the other estimators solutions led to discarding good observations and greater deviations in the estimation.

Small Magnitude Outliers
Starting with the outlier detection, Figure 15 gives the solution classifications for each RE in the 400 scenarios. LS presents its results concentrated in the 'no false positives' classifications ("A", "C" and "E"). This shows a resistance for pointing out false positives in the residuals vector by LS solutions. The IRLS got the highest number of solutions classified as "A", overcoming the other estimators. By checking the results for classifications that identified all outliers ("A" + "B"), the LMS, LTS, and LTS-RC-N2 had the best identification ability. Although the LTA presented many classifications of type "B", this estimator had the lowest number of solutions classified as type "A", similar to SRLS. In Figure 16, we see that LTA presented 1303 false positives in the residuals vector, over the 400 scenarios. It also confirms LS as the estimator with the least amount of false positives (43). In contrast, LS missed most of the blunders, not detecting 572 out of 1000 inserted outliers. Among the REs where IVS was applied, the LMS identified a good part of the outliers without an exaggerated quantity of false positives, missing less than IRLS and LS.

IRLS
LMS LMS-N2  LS  LTA  LTS  LTS-N2  LTS-RC  LTS-RC-N2  SRLS  Identified  641  664  671  428  710  673  657  672  673  674  False  244  553  671  43  1303  616  695  625  615  940  Missed  359  336  329  572  290  327  343  328  327  326  Total  1000  1000  1000  1000  1000  1000  1000  1000  1000  Analyzing the estimator solutions for parameter estimation, Table 2 shows the mean distance of solutions from the official co-ordinates. It also adds the LS result without outliers, to compare the influence of the random errors. Table 2. Mean distance of the solutions for each estimator, organized in the scenarios of 1-4 concurrent small outliers. LS presented the smallest deviations, being the best estimate of the parameters in the case of adopting the solution, without eliminating contaminated observations and making new adjustments. Even in scenarios with four outliers, the LS showed a better result than the REs, although the difference became smaller. The IRLS presented low variation in the solutions, even with an increase of blunders. Whereas most RE presented similar solutions, LTA, in contrast, had the worst estimates.

LS (No Outliers) LS LMS LTS LTS-RC SRLS IRLS LTS-RC-N2 LMS-N2 LTS-N2 LTA
For an application where outliers should be identified, the LS method is not the best. In these cases, according to the results of the experiments, it is recommended to work with IRLS, or, for a better identification with a higher cost in false positives, LMS, LTS, and LTS-RC-N2 obtained more satisfying results.
By analyzing the estimates, LS presented a more solid estimate, even with the lowest detection of blunders. This shows that the outliers which were not detected by LS did not exert great distortion in the estimation. For both scenarios without outliers or schemes that present blunders of small magnitude, LS remains as the best estimator for the parameters. The good redundancy of the network probably contributed to this result. In networks with poorer geometry, the LS will not be as robust to small outliers.

Large Magnitude Outliers
In contrast to the experiments with small outliers, all REs showed greater ease in detecting larger outliers. Most the results were classified as "A" or "B", and rarely as "E" or "F" (see the classifications of the solutions in Figure 17). In this case, LS presented solutions with false positives, "B" and "D" types, resembling the other tested estimators.

LS (No Outliers) LS LMS LTS LTS-RC SRLS IRLS LTS-RC-N2 LMS-N2 LTS-N2 LTA
However, for the final estimate of the co-ordinates, the REs presented more satisfactory estimates. This shows the sensitivity of LS to outliers of great magnitude and the strength of robust methods. Table 4 presents the best estimator, concerning the application and scenario of errors. Table 4. Estimators with best results, according to the demand and scenario. (*) indicates a higher false positives cost.

Conclusions and Future Works
Geodetic networks are the basis for mapping, geoinformation, land cadastre, and other location-based services. They play an important role in society in infrastructural works that depend on accurate control points. This work tested the strategy of using Robust Estimators (REs) in geodetic network adjustment and for detection of outliers.
Several REs were tested. A metaheuristic optimization was conducted for the LTS, LTS-RC, LMS, SRLS, and LTA estimators, whereas, for the IRLS, an iterative process was handled.
The two main contributions of the research were successfully demonstrated: (1) the Independent Vortices Search (IVS) overcame the VS, MVS, HVS, and ABC in all aspects; (2) we performed a deep investigation of several REs, separating the analysis scenarios into no outliers, small outliers, and outliers of great magnitude. This led to the classifications in Table 4, something not explored in the literature before.
In addition, other minor contributions were also presented. One of them was the search space proposed for applying IVS using Equation (27), which is valid not only for geodetic networks, but for any problem based on the Gauss-Markov model. Furthermore, a more detailed analysis of the results, obtained with Equations (7) and (18), was given in Section 5.2, showing which equation is more appropriate in each case.
The experiments with the geodetic network were all performed by applying the IVS, built on the Modified Vortex Search. IVS works with the vortices independently, which can also facilitate the parallelization of procedures which require high performance. This includes exploring more complex problems or even larger geodetic networks.
In the experiments of quality control in the geodetic networks, in situations with small outliers or seeking outlier identification, we do not recommend the LS method. In these cases, according to the experiments, it is better to work with IRLS, or, for a better identification at a cost of more false positives, LMS, LTS, or LTS-RC-N2 can achieve better results. Even though LS detects fewer outliers in these scenarios, it remains as the best estimator for the parameters, as the co-ordinates remain closer to their true values.
For estimating the co-ordinates in scenarios with large-scale outliers, REs present a more satisfactory estimate than LS. This showed the sensitivity of LS to outliers of high magnitude and the strength of the robust methods.
Although the computational cost of REs that use MHs is greater than the classical techniques, it is not possible to establish a cost-benefit at the moment. For this, it is first necessary to consider a minimum amount of FE to get good results, based on the RE and the network dimension. However, a disadvantage of this technique is that is it not possible to estimate the precision of the points. To achieve that, we need an LS estimation after removing the outliers from the observation vector.
Future works can focus on a scalability study and the parallelization of the IVS to achieve better performance. In robust estimation, other strategies can also be studied, such as testing new constraints to REs or checking other REs not contemplated by this work.