A study on multi-GNSS phase-only positioning

Let $\Delta\phi_i\in\mathbb{R}^{f(m-1)}$ denote the vector of double-differenced (DD) observed-minus-computed carrier phase measurements that are collected by two GNSS receivers tracking m common satellites, on f frequencies, at epoch i (i = 1, 2). The distance between the two receivers is assumed short enough so that the DD atmospheric delays can be neglected. With these in mind, the corresponding linearised observation equations read (Teunissen 1997b, Khodabandeh and Teunissen 2018)

$$\begin{equation} \mathsf{E}(\Delta\phi_i) = [\Lambda\otimes I]\,a + [e\otimes G_i]\,\Delta b_i \end{equation} \tag{ 1 }$$

where the unknown DD ambiguities $a\in\mathbb{Z}^{f(m-1)}$ are linked to the measurements through the f × f diagonal matrix $\Lambda = {\textrm{diag}}(\lambda_1,\ldots,\lambda_f)$, with λ_j being the wavelength of the phase data on frequency j (j = 1, ..., f). The 3-vector Δb_i contains the unknown increments of the three-dimensional baseline b_i at epoch i. The corresponding (m − 1)× 3 coefficient matrix is indicated by $G_i = D_m^TA_i$ and is assumed full-rank, where the m × 3 matrix A_i contains the receiver-satellite direction vectors. The m × (m − 1) matrix D_m forms between-satellite differences. The orbital corrections are assumed applied to the observed-minus-computed measurements. In the above equations and in the following, the epoch subscript i is used to emphasise the time dependency of the quantities. In the case of the DD ambiguities a however, the subscript i is omitted. This is because the phase ambiguities are constant in time if no cycle slips occur. It is also important to remark that code data are used to compute the approximate values needed for linearising the observation equations.

Why two epochs: If we stick to only one single epoch of phase data, the system of equations (1) would then be underdetermined as there are f(m − 1) equations in f(m − 1)+3 unknowns, i.e. f(m − 1) ambiguities plus three baseline components. The rank-defect is 3. Assuming that the augmented matrix $[G_1,G_2]$ is of full column rank, at least two epochs of phase data are required to have the model (1) solvable for the unknowns a and Δb_i (i = 1, 2). Under such assumption, the solvability condition of (1) is given by

$$\begin{equation} m \geq 7. \end{equation} \tag{ 2 }$$

Thus phase data of a minimum number of seven satellites is needed. Although phase data of more than two epochs can also be considered, here we confine our study to a two-epoch scenario as it has the potential to realise near real-time phase-only positioning when dealing with high-rate GNSS data, e.g. with 5 Hz or 1 Hz sampling-rates. It should be remarked, for a single-GNSS setup and under unfavourable measurement conditions like urban canyons, that the solvability condition (2) may not always hold (Montenbruck et al 2017).

Why kinematic: In the literature, the phase-only model (1) has been studied under the assumption that the baseline parameters are static over time, i.e. $b_1 = b_2$, see e.g. (Tiberius and de Jonge 1995, Huisman et al 2010, Wang et al 2018). In the present contribution, we relax this stringent assumption and treat the baselines as fully unlinked-in-time parameters. In other words, no restriction is placed on the temporal behaviour of the two baseline vectors b_i (i = 1, 2). The strength of such 'kinematic' model is less than that of its 'static' counterpart in the sense that it has three more unknowns (i.e. the baseline vector of the second epoch). Our goal is therefore to verify whether or not the stated weakness of the kinematic model (1) is of practical relevance for a multi-GNSS landscape where positioning users can track a large number of satellites (∼20), collecting large amounts of phase data on multiple frequencies.

In order to evaluate the variance matrices of the parameter solutions obtained by (1), the stochastic model of the DD phase measurements Δφ_i (i = 1, 2) is assumed given as (Khodabandeh and Teunissen 2015a)

$$\begin{align} Q_{\phi_i \phi_i} & = 2\,\sigma_{\phi}^2\, [I\otimes W_i^{-1}],\quad i = 1,2 \nonumber\\ Q_{\phi_1 \phi_2} & = 0 \end{align} \tag{ 3 }$$

where $W_i = (D_{m}^{T}C D_{m})^{-1}$, with the m × m diagonal co-factor matrix $C = {\textrm{diag}}(w_1^{-1}\!,\ldots,\!w_{m}^{-1})$ whose positive elements w_s (s = 1, ..., m) are the satellite elevation-dependent weights. The zenith-referenced standard deviation of the undifferenced phase measurements is denoted by $\sigma_{\phi}$. The factor '2' in (3) indicates that the variance of the phase measurements is doubled by forming between-receiver differences. According to the second expression of (3), time correlation between the phase measurements Δφ_i (i = 1, 2) is assumed to be absent.

We first present the variance matrix of the parameters of interest, i.e. the baseline parameters, for which we distinguish the following two cases:

• Ambiguity-float: the case that the DD ambiguities are unknown;
• Ambiguity-fixed: the case that the DD ambiguities are successfully resolved as integers.

Lemma.

lemma (Phase-only baseline precision) Let $\hat{b}_i$ and $\check{b}_i$, respectively, be the least-squares ambiguity-float and -fixed solutions of the baseline parameters b_i as given in (1). Their respective variance matrices can be given as

$$\begin{align} Q_{\hat{b}_i \hat{b}_i} & = Q_{\check{b}_i \check{b}_i}\, + \dfrac{2\sigma_{\phi}^2}{f} G_i^{+}\,Q_{\hat{\rho}\hat{\rho}}\,G_i^{+T},\quad i = 1,2 \nonumber\\ Q_{\check{b}_i \check{b}_i} & = \dfrac{2\sigma_{\phi}^2}{f} (G_i^TW_i G_i)^{-1},\quad i = 1,2 \end{align} \tag{ 4 }$$

where $Q_{\hat{\rho}\hat{\rho}} = (\sum\limits_{i = 1}^{2}W_i\mathcal{P}_{G_i}^\bot)^{-1}$, with the projector $\mathcal{P}_{G_i}^\bot = I-G_iG_i^{+}$. The least-squares inverse of G_i is given by $G_i^{+} = (G_i^TW_i G_i)^{-1}G_i^TW_i$.

proof The proof is given in the appendix.

From the expressions (4), one can infer the precision of the baseline solutions obtained by the phase-only model (1). Let us first consider the second expression in which the DD ambiguities are assumed to be fixed to their integers. The variance matrix $Q_{\check{b}_i \check{b}_i}$ is driven by the precision of the phase data $\sigma_{\phi}$, the number of frequencies f and the receiver-satellite geometry through the co-factor matrix $(G_i^TW_i G_i)^{-1}$. The latter is identical to the co-factor of the code-based baseline positioning (Teunissen 1997a), meaning that $Q_{\check{b}_i \check{b}_i}$ is just a downscaled version of its code-based counterpart. Considering that the precision of code-based short-baseline positioning is at the metre-level and that the phase data are almost two orders of magnitude more precise than the code data, the precision of the phase-only ambiguity-fixed solutions $\check{b}_i$ (i = 1, 2) is at the millimetre to centimetre level. Achieving such high-precision solutions relies on the assumption that one successfully maps the unknown ambiguities a to their correct integers, thus requiring an application of ambiguity resolution (Teunissen 1995).

Before applying ambiguity resolution, one should, however, first investigate whether it is worthwhile to do so as the phase-only ambiguity-float solutions $\hat{b}_i$ may already inherit the high-precision of their ambiguity-fixed versions $\check{b}_i$. In the first expression of (4), the ambiguity-float variance matrix $Q_{\hat{b}_i \hat{b}_i}$ is written as the sum of its ambiguity-fixed counterpart $Q_{\check{b}_i \check{b}_i}$ and the positive-definite matrix $(2\sigma_{\phi}^2/f) G_i^{+}\,Q_{\hat{\rho}\hat{\rho}}\,G_i^{+T}$. The latter, in its turn, is driven by the co-factor matrix $Q_{\hat{\rho}\hat{\rho}}$. The inverse of the co-factor matrix $Q_{\hat{\rho}\hat{\rho}}$ is formed as a 'weighted sum' of the projectors $\mathcal{P}_{G_i}^\bot$ (i = 1, 2). This implies, when $G_2\approx G_1$, that the inverse-matrix $Q_{\hat{\rho}\hat{\rho}}^{-1}$ is near rank defect (near singular). This is because the projectors are singular matrices (Teunissen et al 2005). In such cases, the entries of $Q_{\hat{\rho}\hat{\rho}}$ are rather large, indicating that the ambiguity-float solutions $\hat{b}_i$ are poorly estimable. Unfortunately, the approximation $G_2\approx G_1$ holds true for high-rate GNSS data as the receiver-satellite geometries change rather slowly in time. This is due to the high-altitude orbits of the GNSS satellites. The higher the sampling rate, the closer the approximation $G_2\approx G_1$, thereby the poorer the precision of the solutions $\hat{b}_i$ becomes.

To address to what extent the precision of the ambiguity-float baseline solutions is poorer than their fixed versions, we use the concept of gain numbers (Teunissen 1997a). The gain numbers tell us how many times the variance of a function of the baseline solution $\hat{b}_i$ gets smaller after ambiguity-fixing. Since different functions of the baseline can have distinct responses to ambiguity resolution, we employ the average precision gain $\bar{\gamma}$ defined as follows

$$\begin{equation} \bar{\gamma} = \left(\dfrac{|Q_{\hat{b}_i \hat{b}_i}|}{|Q_{\check{b}_i \check{b}_i}|}\right)^{1/3} \end{equation} \tag{ 5 }$$

in which $|\cdot|$ denotes the determinant of a matrix. The following lemma shows the dependency of $\bar{\gamma}$ on the sampling-rate of the phase measurements. Lemma.

lemma (Average precision gain) Let τ be the measurement sampling period, i.e. the reciprocal of the sampling rate. Also let the m × 3 matrix $\dot{A}$ contain the time-derivatives of the entries of $(1/2)\sum_{i = 1}^{2}A_i$ in (1), with $W_2\approx W_1 = W$. Using the first-order approximation $A_2-A_1 \approx \tau\, \dot{A}$ for small values of τ, the average precision gain of the phase-only baseline solutions can be approximated by

$$\begin{align} \bar{\gamma} &\approx \left(\prod\limits_{k = 1}^{3} \left[1+\dfrac{4}{\alpha_k\, \tau^2}\right]\right)^{1/3}\times \dfrac{\left(\prod\limits_{k = 1}^{3} \left[1+\dfrac{\alpha_k\,\tau^2}{4}\right]\right)^{1/3}}{\left(\prod\limits_{k = 1}^{3} [1-\beta_k]\right)^{1/3}}\nonumber\\ &\approx \dfrac{4}{\alpha\,\delta\,\tau^2} \end{align} \tag{ 6 }$$

with $\alpha = \left(\prod\limits_{k = 1}^{3} \alpha_k\right)^{1/3}$ and $\delta = \left(\prod\limits_{k = 1}^{3} [1-\beta_k]\right)^{1/3}$, where the eigenvalues α_k and β_k are, respectively, the roots of the characteristic equations

$$\begin{equation} \begin{array}{ll} \displaystyle |(\dot{G}^TW \dot{G}) - \alpha_k\, (G^TW G)| = 0,&\quad k = 1,2,3\\ \displaystyle |(G^TW \mathcal{P}_{\dot{G}} G) - \beta_k\, (G^TW G)| = 0,&\quad k = 1,2,3 \end{array} \end{equation} \tag{ 7 }$$

with $G = (1/2)\sum_{i = 1}^{2}G_i$, $\dot{G} = D_m^T\dot{A}$, and the projector $\mathcal{P}_{\dot{G}} = \dot{G} (\dot{G}^TW \dot{G})^{-1}\dot{G}^TW$.

proof The proof is given in appendix.

Thus the average precision again $\bar{\gamma}$ is almost inversely proportional to the squared sampling period $\tau^2$ and the positive scalars α and δ. To get some insight into (6), we present time-series of $\sqrt{\alpha}$ and $\sqrt{\delta}$ in figure 1. While the scalar $\sqrt{\alpha}$ is of the order of 10⁻⁴ Hz, the scalar $\sqrt{\delta}$ varies between 0.6 to 1. This implies, for small values of τ, that

$$\begin{equation} \sqrt{\bar{\gamma}}\propto \frac{1}{\tau}\,10^{4} \end{equation} \tag{ 8 }$$

in which τ is in seconds and the symbol $\propto$ means 'is of the order of'. The square-root gain $\sqrt{\bar{\gamma}}$ indicates how many times the standard deviation of the ambiguity-float baseline solutions is roughly larger than its fixed version. Accordingly, when 1 Hz phase measurements are considered (τ = 1 s), the precision of the phase-only ambiguity-float solutions $\hat{b}_i$ is expected to be about tens of metres, i.e. 10⁴ times millimetres. By increasing the sampling period to 10 s (i.e. τ = 10 s), the stated precision improves to several metres which is still way poorer than the sub-centimetre precision of the phase-only fixed solutions $\check{b}_i$.

Zoom In Zoom Out Reset image size

That the phase-only model (1) delivers poorly estimable ambiguity-float baseline solutions may make one inclined to conclude that the model is not applicable to fast precise positioning. This would indeed be the case if successful ambiguity resolution is not realised. The decision whether or not ambiguity resolution is successful is determined by the probability of correct integer estimation, the so-called ambiguity success-rate (Teunissen 1999). The maximum possible success-rate, i.e. the integer least-squares (ILS) success-rate, is governed by the ambiguity variance matrix $Q_{\hat{a}\hat{a}}$. Assuming normally-distributed phase measurements, the ILS success-rate is given by the multivariate integral (ibid)

$$\begin{equation} \textrm{success-rate} = \int_{S_0}\frac{1}{\sqrt{|2\pi Q_{\hat{a}\hat{a}}|}}\exp\left(-\frac{1}{2}||x||_{Q_{\hat{a}\hat{a}}}^2\right)\textrm{d}x \end{equation} \tag{ 9 }$$

in which the ILS pull-in region S₀ is characterised as

$$\begin{equation} S_0 = \{x\in \mathbb{R}^{f(m-1)}|~ ||x||_{Q_{\hat{a}\hat{a}}}^2 \leq ||x-z||_{Q_{\hat{a}\hat{a}}}^2,~\forall z\in\mathbb{Z}^{f(m-1)} \} \end{equation} \tag{ 10 }$$

with the notation $||\cdot||_{Q_{\hat{a}\hat{a}}}^2 = (\cdot)^TQ_{\hat{a}\hat{a}}^{-1}(\cdot)$. According to (9), one only needs the ambiguity variance matrix $Q_{\hat{a}\hat{a}}$ as input to evaluate the ILS success-rate. We therefore present a closed-form expression of $Q_{\hat{a}\hat{a}}$ for the phase-only model (1). Lemma.

lemma (Phase-only ambiguity precision) Let $\hat{a}$ be the least-squares solution of the ambiguity parameters a as given in (1). Then, the ambiguity variance matrix can be given as

$$\begin{align} Q_{\hat{a}\hat{a}} = \sigma_{\phi}^2 [\Lambda^{-1}\mathcal{P}_e^\bot\Lambda^{-1}\otimes W^{-1}]+2\sigma_{\phi}^2 [\Lambda^{-1}\mathcal{P}_e\Lambda^{-1}\otimes Q_{\hat{\rho}\hat{\rho}}] \end{align} \tag{ 11 }$$

with $W \! = \! (1/2)\sum_{i = 1}^{2}\!W_i$ and the projectors $\mathcal{P}_e\! = \!(1/f\,\,) ee^T$ and $\mathcal{P}_e^\bot = I-\mathcal{P}_e$.

proof The proof is given in appendix.

As with the float baseline variance matrices (4), the ambiguity variance matrix $Q_{\hat{a}\hat{a}}$, in (11), also contains the near rank defect co-factor matrix $Q_{\hat{\rho}\hat{\rho}}$. Does this indicate that successful ambiguity resolution is not feasible? To answer this question, we use the Ps-LAMBDA software (Verhagen et al 2013) to numerically evaluate the multivariate integral (9) for a GPS-plus-Galileo scenario. The corresponding ILS success-rates are presented in table 1. The results are shown for different sampling periods τ and numbers of satellites m. The single- and dual-frequency cases refer to the L1(E1) and L1/L5(E1/E5a) signal frequencies, respectively. When a minimum number of satellites (m = 7) is considered, the single-frequency success-rates are almost zero, while the dual-frequency success-rates can be larger than 99% only for the low-rate scenario of τ = 30 s. Upon increasing the number of tracked satellites to m = 20 however, the dual-frequency success-rates can become almost 100% for even high-rate scenarios of τ = 0.2 and 1 s. This highlights the important role played by the number of satellites, and therefore by the multi-GNSS integration, in delivering successful ambiguity-resolved solutions using the phase-only model (1). This notion will be made precise in the following section.

For the sake of comparison, we take the phase-and-code model of single-epoch real-time kinematic (RTK) as reference to analyse the ADOP of the phase-only short-baseline model (1). The corresponding observation equations follow from (1) by sticking to the first epoch i = 1 and including the DD observed-minus-computed code measurements $\Delta p_1\in\mathbb{R}^{f(m-1)}$, that is

$$\begin{equation} \displaystyle\left\{\begin{array}{l}\mathsf{E}(\Delta\phi_1) = [\Lambda\otimes I]\,a + [e\otimes G_1]\,\Delta b_1\\ \displaystyle\mathsf{E}(\Delta p_1) = \quad\quad\quad\,\,\,\,\, +\, [e\otimes G_1]\,\Delta b_1\end{array}\right.. \end{equation} \tag{ 13 }$$

Compare (1) with (13). With the single-epoch model (13), an ambiguity-float baseline solution can be obtained through the second set of the equations, i.e. the code observation equations. The phase-and-code RTK model (13) has been widely used in the literature and its ambiguity resolution performance has been studied for several multi-GNSS scenarios (Odolinski et al 2015, Brack 2017, Paziewski et al 2018). It is therefore important to address how much our phase-only ADOP differs from its rather well-studied RTK counterpart.

Let the positive scalar be the phase-to-code variance ratio. The variance matrix of the DD code measurements Δp₁ follows then by dividing the entries of $Q_{\phi_1 \phi_1}$, in (3), by . The ADOP of the RTK model (13) can be expressed as (Odijk and Teunissen 2008)

$$\begin{align} {\textrm{ADOP}}^{\textrm{RTK}} = \sqrt{2}\,w_o\, \dfrac{\sigma_\phi}{\bar{\lambda}} (1+\frac{1}{\epsilon})^{\frac{3}{2f(m-1)}} \end{align} \tag{ 14 }$$

with $\bar{\lambda} = \prod_{j = 1}^f\lambda_j^{\frac{1}{f}}$ and $w_o = \left(\frac{\sum_{s = 1}^m w_s}{\prod_{s = 1}^m w_s}\right)^{\frac{1}{2(m-1)}}$. To provide insight into (14), let us first make some approximation. The term $\sqrt{2}\,w_o$ is bounded from above by 2 for equal elevation-dependent weights w^s  = 1 (s = 1, ..., m). The precision of current GNSS phase data is also about 1% of their wavelength. Thus $(\sigma_\phi/\bar{\lambda})\approx 0.01$. Since the code data are almost two orders of magnitude less precise than the phase data (i.e.  ≈ 10⁻⁴), the term (1 + 1/) is about 10⁴. We therefore arrive at the following approximation

$$\begin{equation} {\textrm{ADOP}}^{\textrm{RTK}} \approx 0.02\times 100^{\frac{3}{f(m-1)}}. \end{equation} \tag{ 15 }$$

When seven satellites (m = 7) are tracked, the dual-frequency RTK ADOP is about $0.02\times 100^{\frac{1}{4}}\approx 0.06$ cycles which is way smaller than 0.14 cycles. This shows, in the absence of code multipath, that single-GNSS, single-epoch ambiguity resolution is possible with the dual-frequency RTK model (13). By switching to the single-frequency case however, the RTK ADOP is about $0.02\times 100^{\frac{1}{2}}\approx 0.20$ cycles, when m = 7. The single-frequency ADOP reduces to about 0.06 cycles if the number of satellites increases to m = 13. That is why one has to use multi-GNSS data to achieve successful single-epoch ambiguity resolution with the single-frequency RTK model (Odolinski et al 2015).

We now compare the ADOP corresponding to the phase-only model (1) with that of the RTK model (13). The following two phase-only cases are considered:

• Kinematic: the case that the unknown baseline vectors b_i (i = 1, 2) are assumed fully unlinked in time.
• Static: the case that the unknown baseline vector is assumed constant in time, i.e. $b_1 = b_2$;

Lemma.

lemma (Phase-only ADOP compared) Let ${\textrm{ADOP}}^{\textrm{KIN}}$ and ${\textrm{ADOP}}^{\textrm{STA}}$ be the ADOPs corresponding to the kinematic and static phase-only models, respectively. These ADOPs can be linked to ${\textrm{ADOP}}^{\textrm{RTK}}$ in (14) as follows

$$\begin{align} \dfrac{{\textrm{ADOP}}^{\textrm{KIN}}}{{\textrm{ADOP}}^{\textrm{RTK}}} &\approx \dfrac{1}{\sqrt{2}} (\tilde{\epsilon}\,\bar{\gamma})^{\frac{3}{2f(m-1)}} \nonumber\\ \dfrac{{\textrm{ADOP}}^{\textrm{STA}}}{{\textrm{ADOP}}^{\textrm{RTK}}} &\approx \dfrac{1}{\sqrt{2}} \left(\prod\limits_{k = 1}^{3} \tilde{\epsilon}\,\left[1+\dfrac{4}{\alpha_k\, \tau^2}\right]\right)^{\frac{1}{2f(m-1)}}\nonumber\\ \dfrac{{\textrm{ADOP}}^{\textrm{KIN}}}{{\textrm{ADOP}}^{\textrm{STA}}} &\approx (\dfrac{1}{\delta})^{\frac{3}{2f(m-1)}} \end{align} \tag{ 16 }$$

with $\tilde{\epsilon} = \epsilon/(1+\epsilon)$.

proof The proof is given in appendix.

As the first expression of (16) shows, the ADOP and the baseline's average precision gain (5) are closely intertwined. The smaller the average precision gain $\bar{\gamma}$, the smaller the ADOP becomes (Teunissen 1997a). Also note the presence of the factor $1/\sqrt{2}$ in the first two expressions. The factor indicates that the use of two epochs of phase data Δφ_i (i = 1, 2), in the phase-only model (1), can decrease the corresponding ADOP. We now focus our attention on the second expression, i.e. the ADOP of the static phase-only model. With α_k being of the order of 10⁻⁸ Hz² (figure 1), the second expression can be further approximated for small values of τ as

$$\begin{equation} \dfrac{{\textrm{ADOP}}^{\textrm{STA}}}{{\textrm{ADOP}}^{\textrm{RTK}}} \approx \dfrac{1}{\sqrt{2}} \left(\dfrac{2}{\tau}\sqrt{\dfrac{\tilde{\epsilon}}{\alpha}}\,\right)^{\frac{3}{f(m-1)}}. \end{equation} \tag{ 17 }$$

With the approximation $\sqrt{\tilde{\epsilon}/\alpha}\approx 100$, the above ADOP-ratio is about 10 for the single-frequency case, with a minimum number of seven satellites, when τ = 1 s. This ratio reduces to about 2.7 for the dual-frequency case and can be further decreased to only 1.1 when the number of satellite increases to m = 20.

The last expression of (16) enables us to quantify the extent to which the ADOP-performance of the phase-only model (1) improves by switching from kinematic to static mode. Recall, from figure 1, that the scalar $\sqrt{\delta}$ varies between 0.6 and 1. Thus the last expression approaches its maximum value 1, the larger the numbers of frequencies f and satellites m. Figure 2 shows ADOPs corresponding to models (1) and (13) as functions of the number of satellites m for different sampling periods τ. For both the single- and dual-frequency cases, the gap between the ADOPs of the phase-only kinematic and static models gets smaller when more than eight satellites are tracked. The phase-only ADOPs approach their single-epoch RTK versions, the larger the number of satellites. For the low-rate cases (i.e. 10 and 30 s cases), the dual-frequency phase-only ADOPs gets even smaller than their RTK versions when over 15 satellites are tracked. As with the success-rate results in table 1, our ADOP analysis highlights the capability of the dual-epoch phase-only model (1) for delivering high-precision ambiguity-resolved positioning solutions.

Zoom In Zoom Out Reset image size