## 1 Introduction  
Continuous GNSS measurements provide scientists with a time series of coordinate changes not only to detect tectonic signals (lithospheric plate movement, geocentric movement, accumulation of crustal deformation, etc.) but also to respond to temporal variations in the surface load from various sources, such as tidal loading, atmospheric pressure, surface water, etc., as well as systematic errors that have not been fully modeled, including the erroneous modeling of satellite orbits, corrections for the phase center of the satellite/s and receiver antennas, multipath effects, etc. [1]. Both global/regional surface load responses and under-modeled errors lead to the spatial and temporal coherence characteristics of GNSS measurements over a fairly large area. It is known that such coherent characteristics are usually called Common Mode Errors (CMEs), which are the main spatially correlated signals in GNSS measurements [2]. The interaction between surface loading errors and spatial correlation errors can distort the station velocity estimate, increase velocity uncertainty, and even lead to incorrect geophysical interpretation. In order to improve the GNSS signal-to-noise ratio, it is highly recommended to apply various mathematical methods of CME filtering (stacking, reference frame transformation, and statistical signal decomposition techniques) [3, 4, 5, 6]. It should be noted that the results of applying the stackingmethod directly depend on the size of the station networks and the distances between them. The results are satisfactory only if the CME in the regional network is essentially a spatially homogeneous medium. Accordingly, the assumption of equivalent CMEs for all stations in the network, even by introducing appropriate weights, does not solve the problem of identifying stations with strong local effects that may affect the detection of CMEs by this method. A similar situation exists in the case of transforming GNSS solutions into a regional reference system, i.e., stations located at the edge of the region may show distorted signals or some residual CMEs if they are affected by specific local effects. To mitigate the impact of CMEs, statistical signal decomposition techniques are also used. Among them, the most widely used ones are principal component analysis (PCA) [3], multi-channel singular spectrum analysis (MSSA) [5], and independent component analysis (ICA) [4], as well as their numerous modifications and evaluations. Based on the allowable heterogeneous distribution of CMEs (when a CME has only one source) and a more rigorous mathematical structure, these methods have been widely used to remove CMEs from the regional network. The main problem with these methods is the assumption that the residual errors of the GNSS coordinate time series obey a normal distribution law, known as Gaussian law, i.e., the mean and standard deviation are the best estimates used to find the true values and their errors. In the case of a long-term time series from individual GNSS stations, they might contain weak, unaccounted systematic errors [7; 8]. They contribute to the deterioration of the accuracy of the station coordinates, the values of change rates obtained on their basis, as well as the concealment of many weak and transient signals in the coordinate time series (for example, the movement of neighboring faults). The estimation and removal of residual errors are performed using a number of classical statistical analysis methods, such as spatio-temporal filtering, the least squares method (LSM), and maximum likelihood estimation (MLE), which are implemented in all known software packages used for the analysis of the time series of coordinates [9; 10; 11; 12]. In the EPN network, thresholds based on three different categories are used to select the most stable stations: positioning and signal quality, the reliability of determining the rate of change in coordinates, and station stability over the years.