LS-VCE Applied to Stochastic Modeling of GNSS Observation Noise and Process Noise

Abstract

Stochastic models play a crucial role in global navigation satellite systems (GNSS) data processing. Many studies contribute to the stochastic modeling of GNSS observation noise, whereas few studies focus on the stochastic modeling of process noise. This paper proposes a method that is able to jointly estimate the variances of observation noise and process noise. The method is flexible since it is based on the least-squares variance component estimation (LS-VCE), enabling users to estimate the variance components that they are specifically interested in. We apply the proposed method to estimate the variances for the dual-frequency GNSS observation noise and for the process noise of the receiver code bias (RCB). We also investigate the impact of the stochastic model upon parameter estimation, ambiguity resolution, and positioning. The results show that the precision of GNSS observations differs in systems and frequencies. Among the dual-frequency GPS, Galileo, and BDS code observations, the precision of the BDS B3 observations is highest (better than 0.2 m). The precision of the BDS phase observations is better than two millimeters, which is also higher than that of the GPS and Galileo observations. For all three systems, the RCB process noise ranges from 0.5 millimeters to 1 millimeter, with a data sampling rate of 30 s. An improper stochastic model of the observation noise results in an unreliable ambiguity dilution of precision (ADOP) and position dilution of precision (PDOP), thus adversely affecting the assessment of the ambiguity resolution and positioning performance. An inappropriate stochastic model of RCB process noise disturbs the estimation of the receiver clock and the ionosphere delays and is thus harmful for timing and ionosphere retrieval applications.

1. Introduction

Choosing the proper mathematical models, including a functional model, a dynamic model, and a stochastic model, is essential for global navigation satellite systems (GNSS) data processing [1,2]. The functional model describes the relationships between observations and unknown parameters, while the dynamic model specifies the time evolution of unknown parameters [3,4]. A considerable number of studies have established various functional and dynamic models for GNSS positioning, navigation, and timing applications [5,6,7]. However, the stochastic model, which describes the random characteristics of a system, is usually unknown or only partly known. A realistic stochastic model ensures the best linear unbiased estimation, determines an accurate precision description of the unknown parameters, and allows reliable quality control [3,8]. Hence, it is of great importance to establish a realistic stochastic model for GNSS applications.

For stochastic modeling, one must first define the appropriate structure of the covariance matrix. An adequate covariance matrix contains two parts: one for observation noise and another for process noise. Regarding the observation noise, its stochastic properties that reflect the precision of observations, the cross-correlation between different types of observations, and the time correlation of observations can be captured by the covariance matrix [9,10,11]. The stochastic properties of the process noise describe the uncertainty of the dynamic model and should also be properly described by the covariance matrix [8,12]. The correlation between observation noise and process noise is commonly assumed to be zero [8,12].

The covariance matrix is then written as an unknown linear combination of known cofactor matrices [13]. The process of estimating the unknown variance components is referred to as variance component estimation (VCE). Many VCE methods, which tend to differ in estimation principles and distributional assumptions, exist in the literature [14,15,16,17]. This paper focuses on the least-squares variance component estimation (LS-VCE) which is capable of unifying many VCE methods [13,18]. Moreover, it adopts the well-known principle of LS, works with a user-defined weight matrix, and allows one to apply the existing knowledge of LS theories, thus making it simple, flexible, and attractive [13,18].

Many studies have applied the LS-VCE for the stochastic modeling of GNSS observations. Earlier applications of the LS-VCE established a realistic stochastic model for GPS observations [9,19]. Previous studies have also adopted the LS-VCE to analyze the noise characteristics of GPS observations that have been collected by different types of receivers and antennas [10,20,21]. With the emergence of other GNSSs, the LS-VCE contributed to the design of the Chinese BDS a specific stochastic model that considers different variances for geostationary Earth orbit, inclined geosynchronous orbit, and medium Earth orbit satellites [22,23,24]. Galileo multi-frequency observations have also been shown to benefit from a proper stochastic model established by the LS-VCE [25,26,27].

Compared with stochastic modeling of GNSS observation noise, few studies have applied the LS-VCE for the stochastic modeling of process noise. In many applications, knowledge about the process noise for some parameters is still at a rudimentary level. For example, studies in recent years have reported that the receiver code bias (RCB) exhibits short-term variations [28,29,30,31]; however, many studies considered it as a time constant parameter, typically one that is relevant to precise point positioning [32]. RCB is the receiver hardware delay and is composed of contributions from the antenna, cables, amplifiers, and filters. To avoid the adverse effects of RCB variation, one may consider it unlinked in time, with some implementing it in double-differenced real-time kinematic (RTK) positioning, but this results in the model being weaker [33]. The double-differenced RTK positioning eliminates the RCB through differences at each epoch, which is equivalent to modeling the RCB as a time-unlinked parameter. An alternative method is to model the RCB as a random walk process, thus simultaneously ensuring a relatively strong model and capturing RCB variations [34]. However, a realistic stochastic model of the RCB process noise remains to be determined.

Based on the LS-VCE, this paper proposes a method that is able to jointly estimate the variance components for both observation noise and process noise. We first consider the dynamic model as pseudo-observation equations and then integrate them with GNSS observation equations, constructing one covariance matrix that describes both the observation noise and process noise, finally allowing it to be determined through the use of the LS-VCE. To elaborate on the practical aspects of the method, we apply it to jointly estimate the GNSS (GPS+BDS+Galileo) observation noise and the RCB process noise. After the estimates have been validated, we analyze the impact of the stochastic model upon parameter estimation, ambiguity resolution, and positioning.

2. Methods

This section first formulates the general mathematical models that are used for GNSS data processing. The LS-VCE is introduced afterwards. Finally, we present the stochastic modeling method for GNSS observation noise and RCB process noise.

2.1. GNSS Mathematical Models

GNSS mathematical models include a functional model, a dynamic model, and a stochastic model. The functional model describes the relationship between observations and unknown parameters. We write the GNSS functional model as [4]:

$$y\left(i\right)=A_1(i)x_1\left(i\right)+A_2(i)x_2\left(i\right)+\epsilon \left(i\right)$$

(1)

where $y\left(i\right)$ is the observation vector at epoch $i$ ($i=1,2,...,t$), $x_1\left(i\right)$ is the vector for the unknown parameters that can be further modeled by a dynamic model, while vector $x_2\left(i\right)$ contains the unknown parameters that are not able to be modeled by a dynamic model. $A_1(i)$ and $A_2(i)$ constitute the full-rank design matrices, and $\epsilon \left(i\right)$ is the vector of random noise. Noting that, we consider a functional model in which not all of the parameters satisfy a dynamic model; this is a common case in GNSS data processing. For example, a dynamic model is usually unavailable for receiver clocks, while carrier-phase ambiguities are considered as time constants in a continuous arc.

The dynamic model describes the time evolution of the state vector $x_1\left(i\right)$, and it can be written as [4]:

$$x_1\left(i\right)={\mathrm{\Phi }}_{i,i-1}{x}_1\left(i-1\right)+w\left(i\right)$$

(2)

or

$$0=-x_1\left(i\right)+{\mathrm{\Phi }}_{i,i-1}{x}_1\left(i-1\right)+w\left(i\right)$$

(3)

where ${\mathrm{\Phi }}_{i,i-1}$ ($i=2,3,...,t$) is the state transition matrix, and $w\left(i\right)$ is the vector of process noise.

The stochastic model provides the dispersion of random vectors for both the observation noise and process noise, and it can be written as:

$$D\left\{\begin{array}{c}\epsilon \left(i\right)\\ w\left(i\right)\end{array}\right\}=\left[\begin{array}{cc}{Q}_{\epsilon \left(i\right)}& 0\\ 0& Q_{w\left(i\right)}\end{array}\right]$$

(4)

where $D\left\{.\right\}$ denotes the dispersion operator. We describe the dispersion using a covariance matrix that contains the covariance matrix of observation noise $Q_{\epsilon \left(i\right)}$ and the covariance matrix of process noise $Q_{w\left(i\right)}$. The covariance between the observation noise and process noise is assumed to be zero.

When multi-epoch data are collected, we consider the state-transition, Equation (3), to be pseudo-observation equations and combine them with the observation, Equation (1), writing the multi-epoch model as:

$$\begin{array}{l}y=A_1{x}_1+A_2{x}_2+\epsilon \\ 0=B{x}_1+w\end{array}$$

(5)

where $y={\left[y^T\left(1\right),y^T\left(2\right),...,y^T\left(t\right)\right]}^T$ and $A_1=blkdiag\left\{A_1\left(1\right),A_1\left(2\right),...,A_1\left(t\right)\right\}$ in which $blkdiag\left\{.\right\}$ denotes the operator of the block diagonal concatenation of the matrices, $x_1={\left[x_1^T\left(1\right),x_1^T\left(2\right),...,x_1^T\left(t\right)\right]}^T$, $A_2=blkdiag\left\{A_2\left(1\right),A_2\left(2\right),...,A_2\left(t\right)\right\}$, $x_2={\left[x_2^T\left(1\right),x_2^T\left(2\right),...,x_2^T\left(t\right)\right]}^T$, $\epsilon ={\left[{\epsilon }^T\left(1\right),{\epsilon }^T\left(2\right),...,{\epsilon }^T\left(t\right)\right]}^T$, $w={\left[w^T\left(2\right),w^T\left(3\right),...,w^T\left(t\right)\right]}^T$$B={\left[b_2,b_3,...,b_t\right]}^T$, in which $b_i$ ($i=2,3,...,t$) is a $t$ dimensional column vector whose $i_{th}$ element is −1 and whose $i-1_{th}$ element is ${\mathrm{\Phi }}_{i,i-1}$, while the other elements are zero.

We formulate the stochastic model for Equation (5) as:

$$D\left\{\begin{array}{c}\epsilon \\ w\end{array}\right\}=Q_{com}=blkdiag\left\{Q_{\epsilon \left(1\right)},Q_{\epsilon \left(2\right)},...,Q_{\epsilon \left(t\right)},Q_{w\left(2\right)},Q_{w\left(3\right)},...,Q_{w\left(t\right)}\right\}$$

(6)

where $Q_{com}$ is the multi-epoch covariance matrix that captures the stochastic properties for both the observation noise and process noise. Assuming that some elements in the covariance matrix $Q_{com}$ are unknown, we split the matrix $Q_{com}$ into a known part $Q_0$ and an unknown part:

$$Q_{com}=Q_0+{\displaystyle \sum _{k=1}^p{\sigma }_k^2{Q}_k}$$

(7)

where the unknown part contains $p$ variance components ${\sigma }_k^2$ with the deterministic cofactor matrix $Q_k$.

2.2. Least-Squares Variance Component Estimation

Equations (5)–(7) establish the general mathematical models for GNSS data processing. However, part of the stochastic model remains to be determined. Many studies have proposed different variance component estimation methods; here, we focus on the LS-VCE method since it is simple, flexible, and attractive [13,18].

LS-VCE is based on the least-squares principle, and its normal equation can be written as [13]:

$$N{\widehat{\sigma }}^2=L$$

(8)

where ${\widehat{\sigma }}^2={\left[{\sigma }_1^2,{\sigma }_2^2,...,{\sigma }_p^2\right]}^T$ is the vector of the unknown variance components. We determine the entries $n_{kl}$ and $l_k$ of the $p\times p$ matrix $N$ and the $p\times 1$ vector $L$ through:

$$\left\{\begin{array}{c}{n}_{kl}=\frac{1}{2}tr\left(Q_k{Q}_{com}^{-1}{P}_A^{\perp }{Q}_l{Q}_{com}^{-1}{P}_A^{\perp }\right)\\ l_k=\frac{1}{2}{v}^T{Q}_{com}^{-1}{Q}_k{Q}_{com}^{-1}v-\frac{1}{2}tr\left(Q_k{Q}_{com}^{-1}{P}_A^{\perp }{Q}_0{Q}_{com}^{-1}{P}_A^{\perp }\right)\end{array}\right.$$

(9)

where $P_A^{\perp }=I-A{\left(A^T{Q}_{com}^{-1}A\right)}^{-1}{A}^T{Q}_{com}^{-1}$ is an orthogonal projector and $v=P_A^{\perp }{y}_{com}$ is the residual vector. $y_{com}$ and $A$ are written as:

$$y_{com}=\left[\begin{array}{c}y\\ 0\end{array}\right],A=\left[\begin{array}{cc}{A}_1& A_2\\ B& 0\end{array}\right]$$

(10)

which are the observation vector and design matrix in Equation (5), respectively.

Finally, we calculate the unknown variance components through the equation:

$${\widehat{\sigma }}^2=N^{-1}L$$

(11)

In brief summary, Equation (5) combines the functional and dynamic models into a unified framework, making it possible to use one covariance matrix (6) to describe both the observation noise and the process noise. Hence, by means of Equation (7), one can flexibly construct a stochastic model in which the variance components of the observation noise are unknown, those of process noise are unknown, or even both of them are unknown. These unknown variance components can be determined by Equation (11), which is based on the LS-VCE.

2.3. Stochastic Modeling of GNSS Observation Noise and RCB Process Noise

We have already presented the stochastic modeling method of observation noise and process noise. To further elaborate on the method from the practical perspective, we apply the method to estimate the GNSS observation noise and RCB process noise.

2.3.1. Functional and Dynamic Models

To retrieve the stochastic properties of the pure random noise, one has to eliminate the systematic errors. In an ultra-short or zero baseline, the between-receiver single-differenced (SD) model is almost free of systematic errors and is thus suitable for stochastic modeling [24]. Assuming that the between-receiver SD atmospheric delays are eliminated, we write the full-rank single-epoch SD observation equations at frequency $j$ ($j=1,2,...,f$) for $m$ satellites as [35]:

$$\begin{array}{l}{p}_j\left(i\right)=\rho \left(i\right)+e_{m}d\tilde{t}\left(i\right)+e_m{\tilde{d}}_j\left(i\right)+{\epsilon }_{p_j}(i)\\ {\varphi }_j\left(i\right)=\rho \left(i\right)+e_{m}d\tilde{t}\left(i\right)+e_m{\tilde{\delta }}_j\left(i\right)+{\lambda }_j{\tilde{z}}_j\left(i\right)+{\epsilon }_{{\varphi }_j}(i)\end{array}$$

(12)

where $p_j\left(i\right)={\left[p_j^1\left(i\right),p_j^2\left(i\right),...,p_j^m\left(i\right)\right]}^T$ and ${\varphi }_j\left(i\right)={\left[{\varphi }_j^1\left(i\right),{\varphi }_j^2\left(i\right),...,{\varphi }_j^m\left(i\right)\right]}^T$ are the SD code and phase observation vectors, respectively. $\rho \left(i\right)={\left[{\rho }^1\left(i\right),{\rho }^2\left(i\right),...,{\rho }^m\left(i\right)\right]}^T$ is the satellite-to-receiver geometry vector; and ${\tilde{z}}_j\left(i\right)={\left[{\tilde{z}}_j^2\left(i\right),{\tilde{z}}_j^3\left(i\right),...,{\tilde{z}}_j^m\left(i\right)\right]}^T$ is the integer double-differenced ambiguity vector, in which the coefficient ${\lambda }_j$ is the wavelength. $e_m$ is the m-column vector where all of the elements are ones. ${\epsilon }_{p_j}(i)={\left[{\epsilon }_{p_j^1}(i),{\epsilon }_{p_j^2}(i),...,{\epsilon }_{p_j^m}(i)\right]}^T$ is the random noise vector for code observations, while its phase counterpart ${\epsilon }_{{\varphi }_j}(i)$ has the same structure. The parameters denoted by the tildes are not the original parameters but are instead the linear combination of them and are defined in

Table 1. Estimable parameters of between-receiver single-differenced model.

Estimable Parameter	Notation and Interpretation
Receiver clock	$d\tilde{t}\left(i\right)=d{t}_{12}\left(i\right)-d_{12,1}\left(i\right)$
Receiver code bias	${\tilde{d}}_j\left(i\right)=d_{12,j}\left(i\right)-d_{12,1}\left(i\right)\left(j>1\right)$
Receiver phase bias	${\tilde{\delta }}_j\left(i\right)={\delta }_{12,j}\left(i\right)-d_{12,1}\left(i\right)+{\lambda }_j{z}_{12,j}^1\left(i\right)$
Integer ambiguity	${\tilde{z}}_j^m\left(i\right)=z_{12,j}^m\left(i\right)-z_{12,j}^1\left(i\right)$

, where the subscript ‘12’ represents the SD operation between two receivers.

Since the baseline and integer ambiguities can be fixed through the use of long-term observations, we move them to the left side of the equations. Collecting observations of all of the $f$ frequencies yields:

$$y\left(i\right) = \left[\begin{array}{ccc}{e}_{2fm}& I_{2f-1}\otimes e_m& I_{2f}\otimes e_m\end{array}\right]\left[\begin{array}{c}d\tilde{t}\left(i\right)\\ \gamma \left(i\right)\\ \epsilon \left(i\right)\end{array}\right]$$

(13)

where $y\left(i\right)={\left[p^T\left(i\right),{\varphi }^T\left(i\right)\right]}^T$ with the observed-minus-calculated code observation vector $p\left(i\right)={\left[p_1^T\left(i\right),p_2^T\left(i\right),...,p_f^T\left(i\right)\right]}^T$ and its phase counterpart $\varphi \left(i\right)$, $I$ denotes the identity matrix, $\gamma \left(i\right)={\left[{\tilde{d}}_2^T\left(i\right),{\tilde{d}}_3^T\left(i\right),...,{\tilde{d}}_f^T\left(i\right),{\tilde{\delta }}_1^T\left(i\right),{\tilde{\delta }}_2^T\left(i\right),...,{\tilde{\delta }}_f^T\left(i\right)\right]}^T$ is the vector that contains estimable RCBs and receiver phase bias (RPBs), and $\epsilon \left(i\right)={\left[{\epsilon }_{p1}^T(i),{\epsilon }_{p2}^T(i),...,{\epsilon }_{p_f}^T(i),{\epsilon }_{{\varphi }_1}^T(i),{\epsilon }_{{\varphi }_2}^T(i),...,{\epsilon }_{{\varphi }_f}^T(i)\right]}^T$ is the noise vector.

In a multi-epoch case, we integrate the dynamic model with the functional model. We model the estimable receiver clock as white noise and model the receiver code and phase biases as random walk processes, allowing the multi-epoch model to be written as:

$$\begin{array}{l}y=A_1\gamma +A_{2}d\tilde{t}+\epsilon \\ 0=B\gamma +w\end{array}$$

(14)

where we see that Equation (14) has the same structure of Equation (5), but the matrices in Equation (14) are specified for the SD model: $A_1=I_t\otimes e_{2fm}$, $\gamma ={\left[{\gamma }^T\left(1\right),{\gamma }^T\left(2\right),...,{\gamma }^T\left(t\right)\right]}^T$, $A_2=I_t\otimes I_{2f-1}\otimes e_m$, $d\tilde{t}={\left[d{\tilde{t}}^T\left(1\right),d{\tilde{t}}^T\left(2\right),...,d{\tilde{t}}^T\left(t\right)\right]}^T$, and $B={\left[b_2,b_3,...b_t\right]}^T$, in which $b_i$ ($i=2,3,...,t$) is a $t$ dimensional column vector whose $i_{th}$ element is −1 and whose $i-1_{th}$ element is 1, while the other elements are zero.

2.3.2. Formulation of the Stochastic Model

Following the structure in Equation (6), we are able to formulate the covariance matrix for the SD model as:

$$\begin{array}{l}{Q}_{SD}=blkdiag\left\{Q_{\epsilon },Q_w\right\}\\ Q_{\epsilon }=I_t\otimes blkdiag\left\{C_p,C_{\varphi }\right\}\otimes H_m\\ Q_w=I_{t-1}\otimes C_w\end{array}$$

(15)

where $C_p=diag\left\{{\sigma }_{p_1}^2,...,{\sigma }_{p_f}^2\right\}$ and $C_{\varphi }=diag\left\{{\sigma }_{{\varphi }_1}^2,...,{\sigma }_{{\varphi }_f}^2\right\}$ with ${\sigma }_{p_j}^2$ and ${\sigma }_{{\varphi }_j}^2$ are the variances for code and phase observation noise. $C_w=diag\left\{{\sigma }_{d_2}^2,...,{\sigma }_{d_f}^2,{\sigma }_{{\delta }_1}^2,...,{\sigma }_{{\delta }_f}^2\right\}$ with ${\sigma }_{d_j}^2$ and ${\sigma }_{{\delta }_j}^2$ are the variances of the random walk noise for the receiver code and phase biases. $H_m=diag\left\{h_1,...,h_m\right\}$ captures the elevation-dependent precision of the observations, which can be given as [36]:

$$h_m={\left[1+10\exp \left(-\frac{{\theta }_m}{10}\right)\right]}^2$$

(16)

where ${\theta }_m$ is the satellite elevation in degrees.

2.3.3. Estimation of Variances for Observation Noise and Process Noise

To jointly estimate the variances for the observation noise and process noise, we split, following Equation (7), the covariance matrix in Equation (15) into several variance components:

$$Q_{SD}={\displaystyle \sum _{k=1}^{p=4f-1}{\sigma }_k^2{Q}_k}$$

(17)

where ${\sigma }_k^2$ contains $2f$ variance components for the observation noise and $2f-1$ variance components for the process noise of receiver biases. Their cofactor matrices are:

$$\begin{array}{ll}{Q}_k=blkdiag\left\{I_t\otimes blkdiag\left\{0_{m*(k-1)},H_m,0_{m*(2f-k)}\right\},Q_w\right\}& for k=1,...,2f\\ Q_k=blkdiag\left\{Q_{\epsilon },I_t\otimes blkdiag\left\{I_{k-2f-1},1,I_{4f-k-1}\right\}\right\} for& k=2f+1,...,4f-1\end{array}$$

(18)

where $0_{m*(k-1)}$ is a zero matrix of $m$ rows and $\left(k-1\right)$ columns.

After the construction of the normal equations of the LS-VCE by using Equation (9), one can jointly estimate the variances for both the GNSS observation noise and process noise. Experiments will be conducted in the next section, and certain computational problems that may arise will be discussed.

3. Experiments and Results

This section first introduces the experimental setup. We will then analyze the characteristics of the RCBs and RPBs. After showing the estimated variances of the GNSS observation noise and RCB process noise, we validate the estimates and investigate the impact of the stochastic model on parameter estimation in RTK positioning.

3.1. Experiment Setup

Figure 1 shows the experimental setup. The measurement campaign consists of one zero baseline and two ultra-short baselines. Zero (or ultra-short) baselines are commonly used configurations for stochastic modeling since most of the systematic errors (e.g., the multipath) can be eliminated [9,24,37]. In the zero baseline IGG2–IGG3, a Septentrio PolarX5S receiver and a Trimble Alloy receiver were connected to a common antenna. Another antenna was connected to a Septentrio PolarX5S receiver that was 3.2 m away, resulting in the construction of two short baselines: IGG1-IGG2 and IGG1-IGG3. All three receivers collected dual-frequency GPS (L1/L2) [38], Galileo (E1/E5a) [39], and BDS (B1/B3) observations from DOY (day of year) 268 to 270, 2019. The data sampling rate was 30 s. In this experiment, the BDS refers to the third-generation global BeiDou navigation satellite system [40]. We excluded the GLONASS system since its observations are affected by the inter-frequency biases, which should to be calibrated before stochastic modeling [41]. We will use the data on DOY 268 to conduct stochastic modeling and use the data on DOY 269 and DOY 270 to validate the established stochastic models. Studies have demonstrated that one-day data are enough for stochastic modeling [9,24].

3.2. Characteristics of Receiver Biases

Several studies have reported that receiver biases exhibit short-term variations, which are receiver and antenna specific [28,29,30,31]. This section will show the characteristics of the receiver biases that are involved in the data that we collected.

Figure 2 shows the time series of the Galileo RCB. We estimated the RCB in a filtering process and considered the RCB as white noise. The Galileo RCB exhibits variations between –3.6 m to –2.6 m. This confirms the findings that were reported in other studies, which indicated that the RCB varies over the course of a day [28,29,30,31]. Hence, it is unsuitable to model the RCBs as time constants. One may opt to model the RCB as white noise, but this results in the model becoming weaker. One of the methods that is preferable is to model the RCB as a random walk process, provided that the realistic stochastic model of the random-walk noise can be determined.

Figure 3 shows the time series of the GPS L1 RPB. The RPB also shows intra-day variations. Reviewing the definition of RPB that is given in Table 1, the estimable RPB absorbs the RCB at the first frequency. Hence, the RPB time series contains RCB variations and the RCB variations may even dominate the time series. To further investigate this, we turn to Figure 4.

Figure 4 shows the time series of the BDS RCB and the RPB at two different frequencies. The characteristics of the RPB at the two different frequencies are virtually the same. Moreover, the RPB and RCB exhibit similar properties. This further demonstrates that the RCB that is absorbed in the estimable RPB dominates the RPB time series. In other words, all three time series in Figure 4 actually reflect the characteristics of the RCB.

3.3. Variances of GNSS Observation Noise and RCB Process Noise

This section shows the estimated variances of observation noise and RCB process noise. Since the characteristics of the RCB and the RPB at two frequencies were demonstrated to be almost the same, we estimated the same variance for all of them. Since the phase observation precision is much higher than that of the code one, the LS-VCE model may encounter an ill-posed problem. To handle this problem, we followed the strategy in [18] and estimated one same variance for the dual-frequency phase observations. To improve the computational efficiency, we divided the observations into 10-epoch groups, with the final estimates being the means of the 10-epoch estimates [18]. For each group, we estimated the variance component iteratively. The a priori values for code observation noise, phase observation noise, and RCB process noise were set to 0.3 m, 0.003 m, and 0.03 m, respectively. The estimation procedure could be convergent after several (about three to five) iterations.

Figure 5 shows the standard deviations (STDs) of the GPS, Galileo, and BDS SD code observation noise at the first frequency. The STDs of the GPS and Galileo SD code observations are approximately 0.2 m, while the BDS observations are noisier, with STDs of approximately 0.3 m. This indicates that the observation precision differs between systems, which confirms the findings of other studies [10]. The BDS results are absent for a while (UTC hour 12 to 14) since less than four satellites are tracked at that time, making reliable estimates unavailable.

Figure 6 shows the STD of the GPS, Galileo, and BDS SD phase observation noise. The precision of the GPS SD phase observation is approximately 2 millimeters, while the precision of the Galileo SD phase observation can reach 3 millimeters. As for BDS, most of the estimates are under 2 millimeters, indicating that the precision of the BDS phase observations is the highest. The precision of the phase observations is also system specific. Hence, it is critical to establish a realistic stochastic model for each system rather than use one empirical stochastic model for all systems.

Figure 7 shows the estimates of variance components for the RCB process noise. The results indicate that the STD of the RCB process noise is at the sub-millimeter level for all three systems. Such small process noise results in a stronger model compared with a model that considers the RCB as white noise. On the other hand, this random walk noise can capture RCB variations, instead of incorrectly constraining the RCB as a time constant. One may question why the process noise shown in Figure 7 is so small, whereas the RCB time series shown in Figure 2 seems to be much noisier. This is because Figure 7 shows the final estimates with several iterations, while Figure 2 shows the RCB time series calculated by considering the RCB unlinked in time, which are equivalent to the one-iteration results that are not convergent and need further calculations.

Table 2. Estimated STDs of the GPS, Galileo, and BDS observation noise and of the RCB process noise for three baselines.

Baseline	System	RCB (mm)	Code 1 (m)	Code 2 (m)	Phase (mm)
IGG1–IGG2	GPS	0.61	0.21	0.26	2.01
	Galileo	0.69	0.17	0.23	2.57
	BDS	0.74	0.31	0.13	1.73
IGG1–IGG3	GPS	0.62	0.38	0.35	2.35
	Galileo	0.70	0.21	0.26	2.89
	BDS	0.61	0.36	0.17	1.93
IGG2–IGG3	GPS	0.75	0.34	0.27	1.52
	Galileo	0.83	0.14	0.24	2.39
	BDS	0.81	0.31	0.12	1.27

provides all of the estimated STDs of the GPS, Galileo, and BDS observation noise and of the RCB process noise for all three baselines. The results indicate that the precision of observations differs in systems and frequencies. For example, the precision of the Galileo code observations is higher than that of the GPS observations. The BDS B1 code observation is noisier than the BDS B3 code observation. We also see that observations that are collected by different baselines exhibit inhomogeneous precision; this is because the observation quality is related to receivers [20]. The STDs of triple-system RCB process noises are all at the sub-millimeter level at three baselines.

3.4. Validation and Impact of the Stochastic Model

This section aims to validate the stochastic model that was estimated in the last section and to investigate the impact of the stochastic model. We will first focus on the stochastic model of the observation noise and then turn to that of RCB process noise.

To validate the stochastic model of the observation noise, we verify whether the actual precision of RTK positioning matches its formal precision. Here, we refer to the positioning as the ionosphere-fixed RTK in which the between-receiver single-differenced ionosphere delays are assumed to be absent [35]. The actual precision refers to the one calculated by the positioning errors, while the formal precision comes from the covariance matrix that was generated by the least-squares estimator.

Figure 8 shows the RTK positioning using a realistic and an empirical stochastic model. Here, we conducted the RTK positioning in an epoch-by-epoch manner. The empirical stochastic model sets the STDs of code and phase observations as 0.3 m and 3 mm, respectively, and ignores the system-specific and frequency-specific precision. We see that the actual and formal confidence intervals match well when the previously estimated stochastic model is used, indicating that the stochastic model we estimated reflects the realistic precision of observations. However, in the case of an empirical stochastic model, the formal intervals are far away from the actual intervals, indicating that the realistic precision of the observations is higher than the empirical value.

To further investigate the impact of the stochastic model upon ambiguity resolution, Figure 9 shows the ambiguity dilution of precision (ADOP) of the positioning using a realistic and an empirical stochastic model. We see that the ADOP with a realistic stochastic model is much smaller than the ADOP with an empirical one. As the ADOP reveals the formal precision of ambiguities, one can use it as an indicator to decide whether or not to fix the ambiguities or which ambiguities should be fixed. Hence, the unreliable ADOP that is calculated by the unrealistic stochastic model will adversely affect the ambiguity resolution process. We note that, if the empirical values are chosen to be smaller than the realistic values, the ADOP based on the empirical value will be smaller. However, this does not mean that the ambiguity resolution based on the empirical stochastic model will perform better. On the contrary, the smaller ADOP based on the empirical stochastic model indicates that the precision of the ambiguities is overestimated, which would result in incorrect ambiguity resolution [24].

We then validate the STD of RCB process noise that we estimated before. For this purpose, we processed the data in a filtering way and verified whether the residuals of the RCB dynamic model satisfied a zero-mean normal distribution. These RCB residuals were calculated by differencing the RCB estimates in two consecutive epochs. We also verified the consistency between the formal RCB process noise and the actual RCB process noise. The formal RCB process noise refers to the one configured in the filtering process. Here, we set it to 10 m, 0.74 millimeters (estimated before), and 0.001 millimeters, respectively. The actual RCB process noise refers to the STD of the RCB residuals. Figure 10 shows that although the RCB residuals based on three different stochastic models all obey normal distribution, only the results based on our estimated stochastic model exhibit the consistency between formal RCB process noise and actual RCB process noise, which validates the correctness of the stochastic model of the RCB process noise we established.

We then investigate how the stochastic model of the RCB process noise affected the parameter estimation. Figure 11 shows the ADOP, PDOP, and receiver clock precision of the ionosphere-fixed RTK positioning with a realistic and an empirical stochastic model of the RCB process noise. The results indicate that the baseline and ambiguity estimation are insensitive to the RCB stochastic model. However, the RCB stochastic model affects the estimation of the receiver clock, as the formal precision is larger when it is paired with an empirical stochastic model than it is when it is paired with a realistic one. This is because the estimable receiver clock absorbs the RCB, while the baseline and ambiguities are independent of the RCB.

To investigate the impact of the RCB stochastic model upon the ionosphere estimation, we switched the ionosphere-fixed RTK model to the ionosphere-weighted one in which the ionosphere delays are estimated [35]. Figure 12 shows the receiver clock and ionosphere precision in GPS+BDS+Galileo dual-frequency ionosphere-weighted RTK positioning. Except for the receiver clock, the stochastic model of the RCB process noise also affects the ionosphere estimation. Hence, an unrealistic stochastic model of RCB process noise will affect the applications such as timing and ionosphere retrieval.

4. Discussion

In the literature, many studies analyzed the stochastic properties of GNSS observation noise [9,10,20], whereas investigations on the process noise are relatively limited. In this contribution, we proposed a method that is able to jointly estimate the variances of GNSS observation noise and RCB process noise, thereby establishing a realistic stochastic model. By comparing our realistic stochastic model with the empirical stochastic model that is commonly adopted in RTK positioning, we analyzed the impact of the stochastic model on parameter estimation, ambiguity resolution, and positioning.

We focused on the RCB process noise since its stochastic model adopted in many applications is needed to be refined. RCB has long been identified as a time constant parameter in precise point positioning [7,32]. However, many studies that have been published in recent years reported that the RCB exhibits short-term variations [28,29,30]. Our results (in Figure 2, Figure 3 and Figure 4) also confirmed this finding. To avoid the adverse impact of RCB variations, double-differenced RTK positioning considers the RCB as a time-unlinked parameter, but this results in the model becoming weaker [26]. We suggest modeling the RCB as a random walk process that is accompanied by a proper stochastic model. The results (in Figure 7 and Table 2) indicated that the estimated RCB process noise is at the submillimeter level. Applying this realistic process noise into the positioning could create optimal parameter estimation, provided that the stochastic model of GNSS observations is also proper.

Figure 8 and Figure 9 illustrated that the stochastic model of the GNSS observations affects the performance assessment of the positioning and ambiguity resolution. Other studies have demonstrated that an improper stochastic model of observation noise could adversely affect the quality control and deteriorate the positioning performance [24]. As far as the RCB process noise is concerned, previous studies have stated that an improper stochastic model of the RCB process noise could affect parameter estimation [34], while this contribution further clarified which parameters would be affected. As shown in Figure 11 and Figure 12, the stochastic model of the RCB process noise affects the estimation of the ionospheric delay and receiver clock, while the receiver coordinates and ambiguities are insensitive to the stochastic model of RCB process noise.

5. Conclusions

A proper stochastic model is the prerequisite to obtain the best linear unbiased estimation in global navigation satellite systems (GNSS) data processing. This paper proposed a method that is able to jointly estimate the variances for both the observation noise and process noise by using the least-squares variance component estimation (LS-VCE). We applied the proposed method to estimate the observation noise of GPS, BDS, and Galileo, as well as the process noise of the receiver code bias (RCB). We also investigated the impact of the stochastic model upon parameter estimation, ambiguity resolution, and positioning. The main findings and conclusions are discussed below.

The precision of GNSS observations differs in terms of both systems and frequencies. For example, the precision of the Galileo E1 code observation is 0.14 m to 0.21 m, which is higher than that of the GPS and BDS observations at the first frequency. The precision of the BDS B3 code observations is better than 0.2 m, which is higher than its B1 counterpart. The precision of the phase observations also differs, and the BDS phase observations enjoy the highest precision (better than two millimeters) among the three systems. The RCB process noise ranges from 0.5 millimeters to 1 millimeter for all three systems, with a data sampling rate of 30 s.

The stochastic model of observation noise affects the ambiguity resolution and the performance assessment of the positioning. The ambiguity dilution of precision (ADOP) based on the realistic stochastic model is much smaller than that calculated through the use of an empirical stochastic model, indicating that an improper stochastic model affects the predictor for ambiguity resolution performance, which may lead to wrong decisions being made to fix the ambiguities. The formal precision of the positioning matches well with the actual precision when using a realistic stochastic model, while the formal precision is far from the actual precision that is achieved when using an empirical stochastic model. This illustrates that an inappropriate stochastic model affects the assessment of the positioning performance.

The stochastic model of the RCB process noise affects the estimation of the receiver clock and ionosphere delays. The residuals of the RCB dynamic model satisfy a zero-mean normal distribution and the standard deviation (STD) of the residuals is consistent with the process noise that we estimated, thus demonstrating the reliability of the stochastic model. Although the position dilution of precision (PDOP) and ADOP are insensitive to the stochastic model of the process noise, the formal precision of the receiver clock and ionosphere delays are affected. This is because the estimable receiver clock and ionosphere delay absorbs RCB, while the baseline and ambiguities are independent of RCB. Hence, an improper stochastic model of the RCB process noise will adversely affect the timing and ionosphere retrieval applications.

In this paper, we applied the proposed method to jointly estimate the variances of GNSS observation noise and RCB process noise. The method is generally applicable and allows the variances of observation noise and process noise to be estimated according to one’s specific needs and interests. The analysis of this paper is limited to 30 s GNSS observations, one can also use the proposed method to analyze high-rate (e.g., 10 s, 5 s, or 1 s) data.