matrix. Nevertheless, adopting homoscedastic or simplified models of variance and neglecting various correlations between the GNSS signals leads to a reduction in the positioning solution accuracy and reliability. These are crucial issues in relative positioning models using carrier-phase observations, in which correct integer carrier-phase ambiguity resolution (AR) is a necessary condition for precise positioning, since it defines the shape and size of the search space for the ambiguity estimation and contributes to the significance level and power of the validation tests. Unfortunately, nowadays the most common approach is the simplified definition of the stochastic model, which ignores the differences in accuracy and correlation between different GNSS systems. The existing approaches to model the GNSS observation variance components (mostly adapted from the GPS system) use, among others, the following methods: satellite elevation-dependent variance exponential function [1,2] and cosecant variance function [3,4]; Signal-to-Noise Ratio (SNR)-dependent variance model [5, 6, 7], Carrier-to-Noise density Ratio (C/N\({}_{0}\))-based model [8,9] with extensions: SIGMA-e [10,11] and SIGMA-\(\Delta\)[12] models; the ionospheric scintillation-dependent variance model [13]; the receiver tracking error-dependent model (e.g., variance function based on tracking jitters) [13,14]; multipath-dependent variance function [15]; User Ranging Accuracy (URA)-dependent variance function [16,17]. The variance models, as mentioned above, are mainly empirical models based on fitting a function that approximates the variance depending on the selected parameter. Among the methods for determining empirical stochastic models of GNSS observation noise, we can distinguish two main strategies: methods based on modeling observation residues from a zero-length baseline or a short baseline [18, 19, 20] and methods using combinations of both types of measurements [21,22]. An alternative to empirical determination of the VC matrix elements are iterative methods based on the Least-Square Adjustment residuals analysis. These methods use Variance Component Estimation (VCE) algorithms by means of: Minimum Norm Quadratic Unbiased Estimation (MINQUE) [23, 24, 25], Best Invariant Quadratic Unbiased Estimation (BIQUE) [26,27] and Least-Squares Variance Component Estimation (LS-VCE) [28,29]. The iterative methods require an initial VC structure to be assumed, often accomplished using variance functions as mentioned above (e.g., exponential function [30] or cosecant function [4]). Among positioning models using multi-frequency observations or observations from multiple receivers, functional models are most often based on linear combinations of multi-epoch observations (e.g., single or double differences, ionosphere- or geometry-free combinations). In that case, in addition to mathematical correlations [31], the observations physical correlations are also apparent and should be taken into account to accurately capture the variances of observations and the covariances between them in a VC matrix. The physical correlations can be divided into: temporal correlations, mainly due to the atmosphere, multipath and the receiver itself [18,32, 33, 34, 35]; spatial- (between channels) and cross-correlation between observation types (e.g., code and carrier-phase, between frequency) [36, 37, 38]. Determination of correlation parameters, depending on the type of correlation is based on following methods: exponential models [32,39], tropospheric turbulence theory-based models [40,41] and autoregressive moving average (AMRA) modeling of observation residuals [7,42]. Moreover, the VCE algorithms are also applied to study the covariance components of the GNSS observation VC matrix [43, 44, 45, 46]. The models, additionally extended with the characteristics of residual errors (e.g., atmospheric residual errors), constitute a separate group of stochastic models.