## 2 Materials and Methods
### Data Processing and Comparison Method
#### 2.3.3 Data Filtering and Comparisons
This task was executed by implementing the following objective function (Equation (6)) specifically designed to handle the data produced starting from the kinematic survey:  
\[\mathrm{f_{j,k,c}}(\Delta_{\mathrm{D}},\Delta_{\mathrm{H}},\ \
icefrac{{\%}}{{ \mathrm{Q}_{1}}})=\min\left\{\left[\frac{\left(\overline{\Delta_{\mathrm{D}}}{*} {\mathrm{w}_{\Delta_{\mathrm{D}}}}\right)+\left(\overline{\Delta_{\mathrm{H}}} {*}{\mathrm{w}_{\Delta_{\mathrm{H}}}}\right)}{{\mathrm{w}_{\Delta_{\mathrm{D} }}}+{\mathrm{w}_{\Delta_{\mathrm{H}}}}}\right]*\left[\
icefrac{{\%}}{{\mathrm{ Q}_{1}}}\right]^{-1}\right\} \tag{6}\]  
\(\mathrm{f_{j,k,c}}\) is the value of the objective function computed for the k-th (with k potentially ranging from 1 to 9) solution, for the j-th considered satellite combination of four, five, or six satellites, and for the c-th considered constellation (Galileo, GPS, GLONASS); \(\overline{\Delta_{\mathrm{D}}}\) and \(\overline{\Delta_{\mathrm{H}}}\) are the average \(\Delta_{\mathrm{D}}\) and \(\Delta_{\mathrm{H}}\) values for the considered k, j, and c; \({\mathrm{w}_{\Delta_{\mathrm{D}}}}\) and \({\mathrm{w}_{\Delta_{\mathrm{D}}}}\) are the inverse of the variances of the \(\Delta_{\mathrm{D}}\) and \(\Delta_{\mathrm{H}}\) vectors for the considered k, j, and c. Lastly, if \(\
icefrac{{\%}}{{\mathrm{Q}_{1}}}\) was equal to 0, then it was fixed at 0.07%, which was less than one fixed epoch over 1376 kinematic ones. Equation (6) can be considered as the product of a weighted mean of the planimetric and altimetric deviations multiplied by the inverse of the percentage of fixed solutions. It was designed to create a trade-off among solutions, considering the goodness of the deviations (thus the discrepancies from the reference trajectory) and the quality of the solutions (fixed or float). In this way, if a computed GS14 trajectory was characterized by good deviations towards the reference trajectory (thus, small \(\overline{\Delta_{\mathrm{D}}}\) and \(\overline{\Delta_{\mathrm{H}}}\)) but with a low percentage of fixed solutions, there were few probabilities that it would be chosen when compared to solutions with slightly different deviations but characterized by a higher fix percentage. After the application of the objective function (Equation (6)), 48 vectors for each of the 1376 kinematic epochs, were created to perform the comparisons:  
1. 16 vectors for the 4, 5, and 6 satellites combinations;
2. Among the 16 vectors for the considered combination (e.g., of 4 satellites), 8 vectors were generated for Galileo and GPS comparisons, and 8 vectors were generated for Galileo and GLONASS comparisons;
3. Considering, for example, Galileo and GPS comparisons (the same is valid for Galileo and GLONASS comparisons), the 8 vectors were produced by 4 couples of vectors finally used for the comparisons;
4. The 4 couples were built using, for each considered epoch (t), all the values found for the considered epoch t for planimetric and altimetric deviations, separated for fixed and float solutions. That is (e.g., for Galileo and GPS comparisons): * (I): two vectors \(\Delta_{\mathrm{D,t}}\) referred to the t-th epoch in case of fixed solutions (one for Galileo and one for GPS system); * (II): two vectors \(\Delta_{\mathrm{D,t}}\) referred to the t-th epoch in case of float solutions (one for Galileo and one for GPS system); * (III): two vectors \(\Delta_{\mathrm{H,t}}\) referred to the t-th epoch in case of fixed solutions (one for Galileo and one for GPS system); * (IV): two vectors \(\Delta_{\mathrm{H,t}}\) referred to the t-th epoch in case of float solutions (one for Galileo and one for GPS system);
5. Particularly, for each considered couple of vectors (e.g., couple I) related to a specific epoch t, the mean values for the Galileo analyzed parameter and for the same GPS (or GLONASS) parameter were computed and compared;
6.