# Uncertainty Quantification of First Fix in a Time-Differenced Carrier Phase Observation Model

Hakim Cherfi, Julien Lesouple, Joan Solà, Paul Thevenon

► To cite this version:

Hakim Cherfi, Julien Lesouple, Joan Solà, Paul Thevenon. Uncertainty Quantification of First Fix in a Time-Differenced Carrier Phase Observation Model. Sensors, 2025, 25 (11), 10.3390/s25113480. hal-05096928

HAL Id: hal-05096928

https://hal.science/hal-05096928v1

Submitted on 4 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Uncertainty Quantification of First Fix in a Time-Differenced Carrier Phase Observation Model

Hakim Cherfi <fnref n="1" />,* , Julien Lesouple <fnref n="1" />, Joan Solà <fnref n="2" /> and Paul Thevenon <fnref n="1" /> ID

**Abstract:** This paper presents an uncertainty quantification analysis of the first fix in a time-differenced carrier phase (TDCP) observation model. TDCP is a widely used method in GNSS-based odometry for precise positioning and displacement estimation. A key point in the TDCP modeling is the assumption that the GNSS receiver's initial position is perfectly known, which is never exactly the case in real-world applications. This study assesses the impact of initial position errors on estimated displacement by formulating a correct TDCP model and a misspecified one, where the first position is not correct. Theoretical derivations provide a generic framework of estimation under the misspecified model and its associated mean squared error (MSE), as well as estimation performance bounds through the misspecified Cramer Rao bound (MCRB) for the considered case. These theoretical considerations are then used to build an estimator of the receiver's displacement, with comparisons to the MCRB for performance evaluation. Extensive simulations using realistic GNSS geometry assess the influence of a first-fix error under various conditions, including different time intervals, first-fix error norms, and first-fix error direction. As an example, it is shown that for the considered geometry, if a TDCP of $t_2 - t_1 = 1$ s is built with an initial first fix error norm $||\Delta r_1|| = 10$ m, then it introduces an estimation of the displacement, with an error of norm equal to 1.3 mm, at most. The results indicate that the displacement estimation error is linearly related to the initial position error and the time interval between observations, highlighting the importance of accurate first-fix estimation for reliable TDCP-based odometry. The findings contribute to highlighting the order of magnitude of errors on solutions as a function of the error on parameters.

**Keywords:** GNSS; TDCP; uncertainty quantification; odometry; positioning error; navigation

<figure><img src="image_2.png" /></figure>

### check for updates

Academic Editor: Baocheng Zhang

Received: 16 April 2025

Revised: 16 May 2025

Accepted: 30 May 2025

Published: 31 May 2025

**Citation:** Cherfi, H.; Lesouple, J.; Solà, J.; Thevenon, P. Uncertainty Quantification of First Fix in a Time-Differenced Carrier Phase Observation Model. *Sensors* **2025**, *25*, 3480. https://doi.org/10.3390/s25113480

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Odometry is a fundamental technique used in robotics, autonomous vehicles, and geospatial applications for estimating the change in position of a moving system. It relies on sensors or data from systems like wheel encoders, inertial measurement units, vision sensors, lasers, or global navigation satellite systems (GNSS) to measure the motion of a vehicle or object. GNSS-based odometry has become an essential tool for achieving precise and reliable positioning and displacement estimation over long distances.

In GNSS-based odometry, the system estimates its positions and/or displacements by processing satellite observations. One of the most precise observations is the carrier phase observation. The carrier phase refers to the phase of the carrier wave used by satellites to transmit signals. This carrier phase observation is particularly valuable because it provides

<footnote n="1">¹ Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, 7, Avenue Edouard Belin, 31400 Toulouse, France; julien.lesouple@enac.fr (J.L.); thevenon@recherche.enac.fr (P.T.)</footnote>

<footnote n="2">² Institut de Robòtica i Informàtica Industrial (IRII), CSIC-Universitat Politècnica de Catalunya, 08028 Barcelona, Spain; jsola@iri.upc.edu</footnote>

<footnote>* Correspondence: hakim.cherfi@enac.fr</footnote>
much higher accuracy (on the order of millimeters) compared to standard pseudorange
measurements (with meter-level accuracy), making it ideal for precise positioning and
odometry. However, this observation is ambiguous, the ambiguity being the integer number
of wavelengths traveled by the signal. One way to get rid of this ambiguity, if constant over
time, is to use the subtraction of two consecutive carrier phase observations and form what
is known as the time-differenced carrier phase (TDCP) observation. This method is widely
used in GNSS-based odometry to estimate the change in position over time. In TDCP, it is
assumed that the integer ambiguities present in both carrier phase observations cancel out
when the subtraction is built, allowing for highly accurate displacement estimation.

Uncertainty quantification refers to the process of quantifying the uncertainties in a
quantity of interest by propagating the uncertainties in the input of a system. Uncertainty
in GNSS measurements arises from various factors, such as observation noise or an approx-
imate assumed position that is used when building observation models. Reviews about
sources of error and their impact within GNSS systems are available in [1].

This paper proposes to focus on the TDCP model: when constructing the TDCP
equations as a function of the receiver displacement, the receiver's position must be known
at one of the two epochs. This knowledge on position is never perfect and can affect the
estimated displacement. To account for this, uncertainty quantification is used to assess
the impact of this inaccuracy on the displacement estimate. By quantifying the uncertainty
in the receiver's position, one can evaluate how this error propagates through the TDCP
equations and influences the displacement estimation. This helps in understanding the
potential error margins and ensuring that the odometry results remain reliable, even when
position uncertainties are present, thus improving the reliability of GNSS-based odometry
in real-world applications.

## 2. Related Work

### 2.1. TDCP Models

Using a GNSS receiver, a system could estimate its speed vector at a given instant by
means of Doppler observations. Once integrated, these speeds provide an estimation of the
system's displacement. It is known that the Doppler observations, when integrated, are noisier
than carrier phase observations, though not sensitive to cycle slip [2]. For this reason, it is
interesting to estimate a receiver displacement between two instants by means of TDCP.

In [3], the authors assume a perfectly known receiver position at the first epoch and
use this information to run a Kalman filter with TDCP and substracted TDCP observations
between satellites, along with a cycle slip detector. The source of the first position is not
specified in the experiments. Also, the receiver-to-satellite unit vectors at the two epochs
are known. To assess potential cycle slips, the cycle slip detector first generates a set of valid
solutions (from a previous estimate, integrated IMU measurements and IMU specifications),
and checks if the TDCP solution is within that set.

In [3,4], the authors build substracted TDCPs, i.e., the subtraction of TDCPs between
two different satellites, in order to remove the receiver clock bias difference. The subtraction
is done with the satellite with the highest elevation, as the receiver most likely receives a
direct signal from this one. Still, to estimate the receiver displacement with these subtracted
TDCPs, the authors assume known receiver-to-satellite vectors and receiver position at the
first epoch of substracted TDCPs.

Some authors go further in the approximations of the TDCP model. In [5,6], the
receiver-to-satellite unit vectors are assumed constant at the two TDCP epochs, leading to a
unique receiver-to-satellite vector. It should be noted that this approximation requires the
knowledge of the receiver position. Then, the receiver and satellite displacement vectors
are projected onto this vector thanks to a scalar product. This makes the TDCP a linear model of the receiver displacement.

The same approximation is made in [2,7], where the author explains that the TDCP model is linearized around the estimated position at the previous epoch. Also, the authors build the single difference of carrier phase observations between satellites and between epochs.

In [8], a similar TDCP model is built, with a time between the two epochs ranging from 1 s to 60 s. Regarding the geometric part, the model includes two terms, the first one being a change of range between the satellite and the receiver, which requires the receiver position, and a second term where the estimated receiver displacement vector is projected onto the unit line-of-sight vector at the second epoch, which also requires the receiver position to be known. Also, a single line-of-sight vector is considered, even if the time of interval is 60 s. This TDCP model is used to build a TDCP factor inside a factor graph optimization problem. At the first epoch, a first position is set from a prior factor, and for other epochs, the first position is set from previous nodes. It should be noted that cycle slip factors are used to handle cycle slips.

In [9], the geometric part of the TDCP model is split into three terms, the first one due to the satellite motion along the line-of-sight, a second term that contains the receiver position at the first epoch, and the last term is a linear term of the receiver displacement and of the line-of-sight vector at the second epoch. All these terms require the knowledge of the receiver coordinates at the two epochs.

In [10], the authors leave the TDCP model as a function of the two receiver positions. The TDCPs are used in a factor graph optimization problem, where relative displacement factors from preintegrated inertial measurements are used. In the experimental part of their work, their proposed method outperforms two filtering methods that also use TDCP observations. This is likely due to a prior factor that sets the receiver position at the first epoch, although no details are given regarding this assumption.

Overall, within the GNSS community, many approximations of the TDCP model are provided, which all assume the known receiver coordinates in order to estimate a displacement.

## 2.2. Uncertainty Quantification

Uncertainty quantification can be defined as the process of assessing (and if possible, reducing) the impact of an uncertainty within a system on its output. It is a very general topic that applies to any engineering problem. In the case of a deterministic system, the uncertainty can either come from a lack of knowledge or from a probabilistic modeling of it. Indeed, when the system is probabilistic, the uncertainty directly comes from it. Regarding the quantification, multiple methods exist, depending on where the uncertainty comes from.

In [11], the authors deliver a Python library for uncertainty quantification, for very general cases, ranging from physics to applied mathematics, and with many algorithms.

Within the GNSS community, uncertainty quantification has already been applied. For example, in [12], the authors consider multiple sources of uncertainty, which propagates in GNSS reflectometry observations to estimate sea level. The quantification is made using the RMSE metric. In [13], the authors estimate a GNSS receiver position using noisy observations. The uncertainty is modeled by an Ornstein-Uhlenbeck random process, a Brownian random process, with learned parameters. The observations are run in a Kalman filter, and they find out that all these methods outperform, in terms of uncertainty, the direct approach which consists of running the Kalman Filter in a traditional way. In [14], from confidence volumes of Gaussian distributions, the authors provide analytical extensions of these volumes into zonotopes that take into account uncertainties of an observation model containing potential biases and/or random errors.
## 2.3. Overview of Proposed Method

In the estimation community, work on estimators built with misspecified models exist [15–17]. On the other hand, in the GNSS community, to the best of the author’s knowledge, the impact of a misspecified first position in a TDCP model has never been studied. Therefore, a correct TDCP model, as well as a misspecified one, is going to be built. The misspecified one will be different from the correct one by the fact that the true receiver position at the first epoch is never known perfectly. The receiver displacement will be estimated, and the receiver position at the second epoch will be deduced from the receiver position at the first epoch and the estimated displacement. No approximations will be made in these models. They will contain the subtraction of ranges, but these ranges will not be projected onto a line-of-sight vector, as has been done by previously mentioned authors. These two models will be compared experimentally, with multiple values around the true receiver position at the first epoch, in order to assess the importance of the receiver position at the first epoch in a TDCP model to properly estimate a receiver displacement.

# 3. Carrier Phase Observation Model and TDCP

## 3.1. Carrier Phase Observation Model

A GNSS receiver at epoch $t$ of coordinates $\mathbf{r}_t \in \mathbb{R}^3$ with a clock bias with respect to GPS time $dt_t^r$ is tracking a satellite of coordinates $\mathbf{s}_t \in \mathbb{R}^3$ with a clock bias with respect to GPS time $dt_t^s$.

The carrier phase observation is modeled as follows:

$$ \lambda\Phi_t = \|\mathbf{s}_t - \mathbf{r}_t\| + c \cdot (dt_t^r - dt_t^s) + \lambda N_t - I_t + T_t + m_t + \varepsilon_t^\Phi \quad (1) $$

where the subscript $t$ is used to indicate that the parameters are expressed at epoch $t$, $\Phi_t$ is the carrier phase raw measurement, $\lambda$ is the wavelength of the considered frequency band, $c$ is the speed of light, $I_t$ is the ionospheric delay, $T_t$ is the tropospheric delay, $m_t$ is the multipath error, $N_t \in \mathbb{Z}$ is the integer ambiguity, $\varepsilon_t^\Phi$ is a centered Gaussian noise with standard deviation $\sigma_\Phi$, and the norm of a vector $\mathbf{v} \in \mathbb{R}^3$ is defined as $\|\mathbf{v}\| = \sqrt{\mathbf{v}^T \mathbf{v}}$.

In this model, the satellite coordinates $\mathbf{s}_t$ and clock bias $dt_t^s$ are known from the ephemeride information. In a precise point positioning (PPP) framework, the unknown quantities that are estimated are the receiver coordinates $\mathbf{r}_t$, the clock bias $dt_t^r$, the integer ambiguity $N_t$, and possible other atmospheric states.

## 3.2. TDCP Observation Model and Corrections

Considering two carrier phase observations as in Equation (1) at epochs $t_1$ and $t_2$ from the same satellite, the TDCP model is built:

$$ \lambda(\Phi_2 - \Phi_1) = \|\mathbf{s}_2 - \mathbf{r}_2\| - \|\mathbf{s}_1 - \mathbf{r}_1\| + c \cdot (dt_2^r - dt_2^s) - c \cdot (dt_1^r - dt_1^s) \\ + \lambda \cdot (N_2 - N_1) - I_2 + I_1 + T_2 - T_1 + m_2 - m_1 + \varepsilon_2^\Phi - \varepsilon_1^\Phi. \quad (2) $$

Correcting all the atmospheric terms [18,19] and satellite clock biases, assuming the difference between the two multipath terms is an additive centered Gaussian noise, and assuming no cycle slip ($N_2 - N_1 = 0$), the model (2) becomes:

$$ \lambda(\tilde{\Phi}_2 - \tilde{\Phi}_1) = \|\mathbf{s}_2 - \mathbf{r}_2\| - \|\mathbf{s}_1 - \mathbf{r}_1\| + c \cdot (dt_2^r - dt_1^r) + \varepsilon \quad (3) $$

where $\tilde{\Phi}_t$ corresponds to the corrected carrier phase observation, $\varepsilon$ is the contribution of the two independent errors $\varepsilon_2^\Phi - \varepsilon_1^\Phi$, as well as all other mutually independent Gaussian centered residual errors after correction. Therefore, $\varepsilon$ is a centered Gaussian noise with
standard deviation $\sigma_{TDCP} > \sqrt{2}\sigma_{\phi}$. The assumption of constant ambiguity requires the use of a properly working cycle slip detector [20].

Written this way, the model in Equation (3) is a function of the receiver coordinates $r_1$ and $r_2$ and receiver clock biases $dt_1^r$ and $dt_2^r$. While the use of multiple observations from other satellites (at least 4 for each epoch) with this model would still permit one to estimate all of these quantities, depending on the time elapsed between epochs $t_1$ and $t_2$, the obtained result would be highly sensitive to noise, and the mean error of an optimal unbiased estimator would be very large. This is due to the fact that between two epochs, the receiver-to-satellite unit vectors would not change much, hence a linear approximation that is badly conditioned. This point discourages a user from estimating all these quantities directly using TDCP observations unless considering a very large time interval between the two considered epochs.

Regarding the ephemerides, in [9], the authors compared the use of final products (i.e., high-precision satellite coordinates in post-processing conditions) versus real-time navigation data and figured out that the use of final products would not improve the solution with respect to the use of navigation data, as between the two epochs, if the same set of ephemerides is used, the satellites positions errors are approximately constant. On the other hand, they figured out that when the set of ephemeride parameters is updated, a gap appears in the satellite positions, which adds an important error in a TDCP model if this update happens between the two epochs. To cancel this phenomenon, it is therefore important to use the same set of ephemerides in a TDCP model.

### 3.3. TDCP Model with True Receiver Position and Receiver Displacement

The model in Equation (3) can be expressed in terms of receiver displacement and clock bias difference. The receiver displacement vector between epochs $t_1$ and $t_2$ is defined as $r_{1,2} \triangleq r_2 - r_1$, and the receiver clock bias difference between epochs $t_1$ and $t_2$ is defined as $dt_{1,2}^r \triangleq dt_{t_2}^r - dt_{t_1}^r$. With these notations, the model in Equation (3) is written as:

$$ \lambda(\tilde{\Phi}_2 - \tilde{\Phi}_1) = f_{r_1}(r_{1,2}, dt_{1,2}^r) + \varepsilon. \quad (4) $$

where

$$ f_{r_1}: \mathbb{R}^3 \times \mathbb{R} \to \mathbb{R} \\ (x_1, x_2) \mapsto \|s_2 - (r_1 + x_1)\| - \|s_1 - r_1\| + c \cdot x_2. \quad (5) $$

The index indicates that this function is not a function of the receiver coordinates but of its displacement, with the position at epoch $t_1$ set to $r_1$, and not of a receiver clock bias, but of its clock bias difference. Given a set of visible satellites, this observation model and TDCP observations can be built for each satellite. The set of models and observations can then be used to estimate the receiver displacement vector and clock bias difference by solving the maximum likelihood problem. In practice, the receiver position $r_1$ is not perfectly known. Therefore, in the next subsection, an error term will be introduced, leading to a different model.

### 3.4. TDCP Model with Assumed Receiver Position and Receiver Displacement

The previous model $f_{r_1}$ in Equation (5) assumed perfect knowledge of the true receiver coordinates $r_1$. In practice, these coordinates can be inaccurate. A first fix denotes an estimated initial position $\hat{r}_1$ of the true initial position $r_1$. Therefore, the error between the true coordinates and the assumed ones are defined as $\Delta r_1 \triangleq \hat{r}_1 - r_1$. From these notations, the true model is the one in Equation (5) while the assumed one is:
$$f_{\hat{r}_1}: \mathbb{R}^3 \times \mathbb{R} \to \mathbb{R} \\ (x_1, x_2) \mapsto \|s_2 - (\hat{r}_1 + x_1)\| - \|s_1 - \hat{r}_1\| + c.x_2. \qquad (6)$$

One can notice that $\hat{r}_1 = r_1 \Leftrightarrow f_{\hat{r}_1} = f_{r_1}$.

In order to quantify the difference between the two models as a function of $\Delta r_1$, one defines the difference between the assumed model Equation (6) and the true model Equation (5), expressed as a function of the error $\Delta r_1$, for the true receiver displacement and the true receiver clock bias difference:

$$\Delta f: \mathbb{R}^3 \to \mathbb{R} \\ \Delta r_1 \mapsto f_{r_1+\Delta r_1}(r_{1,2}, dt_{1,2}^r) - f_{r_1}(r_{1,2}, dt_{1,2}^r). \qquad (7)$$

It can be noted that this application cancels all the clock biases terms, and only the difference in the geometric ranges remains:

$$\Delta f(\Delta r_1) = \|s_2 - (r_1 + \Delta r_1 + r_{1,2})\| - \|s_1 - (r_1 + \Delta r_1)\| \\ - \|s_2 - (r_1 + r_{1,2})\| + \|s_1 - r_1\|. \qquad (8)$$

Therefore, the first fix $\hat{r}_1$ may have an impact on the estimated receiver displacement and clock bias difference.

In general, the use of a misspecified model to estimate parameters has an impact on the estimation. The objective of the following sections is to experimentally assess the estimation error due to a fix error $\Delta r_1$.

### 3.5. Estimation of a Receiver Displacement and Clock Bias Difference from TDCP Observations

Now that a misspecified TDCP model has been written in Equation (6), one can estimate a receiver displacement and clock bias difference from a set of $m$ visible satellites at the two epochs $t_1$ and $t_2$ and from the set of TDCP observations at these epochs. A receiver's displacement estimator is going to be defined and is denoted $\hat{r}_{1,2}$. The second epoch receiver's position estimator is denoted $\hat{r}_2 \triangleq \hat{r}_1 + \hat{r}_{1,2}$. A 2D representation of the TDCP model, for a single satellite, is visible in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1. Representation of the TDCP model in 2D. The blue disks represent the satellite positions at the two epochs, the green crosses represent the true receiver positions, the red ones represent the assumed receiver positions.</figcaption></figure>
To assess the degradation of the displacement estimation due to a misspecified first fix, one can explicitly write this estimator, compute its statistical properties, and assess the impact of Δr₁ on it. In the next section, theoretical elements are provided to achieve these objectives.

## 4. Estimation: Computation of Mean Squared Error and Its Lower Bound

This section uses a more general framework than the one used so far in this article and provides the theoretical elements to build an estimator of a vector of parameters living in a Euclidean space from a multivariate Gaussian random variable. In the considered framework, the mean value of the Gaussian random variable is misspecified. This uncertainty is going to impact the estimator's performances, and one proposes to quantify this impact. Additionally, this section provides a theory to compute the statistical properties of this estimator. Finally, the link with the original problem (the estimation of a GNSS receiver's displacement and clock bias difference, from a set of TDCP observations) is made.

### 4.1. Estimator of a Real Vector

Given an unknown vector $\theta \in \mathbb{R}^n$ and a function $g : \mathbb{R}^n \to \mathbb{R}^m$ differentiable at $\theta$, one obtains a vector of observations:

$$y = g(\theta) + \varepsilon, \varepsilon \sim \mathcal{N}(0, R). \qquad (9)$$

where $\mathcal{N}(\mu, \Sigma)$ denotes the multivariate Gaussian distribution with mean vector $\mu$ and covariance matrix $\Sigma$. In this study, the function $g$ is unknown, i.e., the construction of an estimator of $\theta$ will be based on a function $h : \mathbb{R}^n \to \mathbb{R}^m$.

Since $y \sim \mathcal{N}(g(\theta), R)$, but $g$ is unknown, one defines the misspecified likelihood function:

$$q_y : \mathbb{R}^n \to \mathbb{R} \\ x \mapsto \frac{1}{(2\pi)^{\frac{m}{2}} \sqrt{|R|}} \exp\left(-\frac{1}{2} \|y - h(x)\|_R^2\right) \qquad (10)$$

where $\|v\|_R^2 \triangleq v^T R^{-1} v$. The mismatched maximum likelihood estimator (MMLE) is defined as:

$$\hat{\theta} \triangleq \arg \max_{x \in \mathbb{R}^n} q_y(x) = \arg \min_{x \in \mathbb{R}^n} \|y - h(x)\|_R^2 \qquad (11)$$

as long as the above equation is a measurable function of $y$. $\hat{\theta}$ is Gaussian if and only if $h$ is a linear or affine function. In general, computing its distribution or statistical properties is not straightforward. Instead, an approximation of $\hat{\theta}$ is defined: assuming that the function $h$ is differentiable on a set $S \subseteq \mathbb{R}^n$, then one can choose a vector $\theta_0 \in S$, and the function can be written:

$$\forall x \in \mathbb{R}^n, h(x) = h(\theta_0) + H_{\theta_0}(x - \theta_0) + o(\|x - \theta_0\|) \qquad (12)$$

where $H_{\theta_0}$ represents the Jacobian matrix of $h$ at $\theta_0$. This way, the function $h$ is replaced by its first order development at $\theta_0$ in the misspecified likelihood Equation (10), leading to:

$$\hat{\theta} \approx \arg \min_{x \in \mathbb{R}^n} \|y - (h(\theta_0) + H_{\theta_0}(x - \theta_0))\|_R^2. \qquad (13)$$

The vector $\theta_0$ can be obtained by running the Gauss-Newton algorithm with a sample of $y$ and an arbitrary initial estimate. Under certain conditions on the function $h$ and on the
step size [21], the Gauss-Newton algorithm is guaranteed to converge to a local minima and give an estimate (denoted $\theta_0$). The proposed approach assumes that the function $h$ is not too nonlinear (in particular, that from the definition of an operator norm, $|||H_{\theta_0} - H_{\theta}||| \approx 0$).

If $H_{\theta_0}$ represents an injective application, one defines the matrix:

$$S \triangleq (H_{\theta_0}^T R^{-1} H_{\theta_0})^{-1} H_{\theta_0}^T R^{-1}. \qquad (14)$$

The solution of the optimization problem Equation (13) is given by:

$$\hat{\theta} \approx S(y - h(\theta_0)) + \theta_0. \qquad (15)$$

In this case, $\hat{\theta}$ is approximated by a Gaussian random vector.

## 4.2. Properties of the Estimator

Now that an estimator Equation (15) is available, some of its statistical properties are provided. One determines that:

* Since $y \sim N(g(\theta), R)$, the estimator's mean vector is $E[\hat{\theta}] = S(g(\theta) - h(\theta_0)) + \theta_0$;

* The estimator's bias is $b[\hat{\theta}] = S(g(\theta) - h(\theta_0)) + \theta_0 - \theta$;

* Since $\hat{\theta} - E[\hat{\theta}] = S\epsilon$, the estimator's covariance matrix is $V[\hat{\theta}] = SRS^T = (H_{\theta_0}^T R^{-1} H_{\theta_0})^{-1}$.

Defining the mean squared error operator $MSE[\hat{\theta}] \triangleq E[(\hat{\theta} - \theta)(\hat{\theta} - \theta)^T]$, it can be shown that:

$$MSE[\hat{\theta}] = V[\hat{\theta}] + b[\hat{\theta}]b[\hat{\theta}]^T. \qquad (16)$$

## 4.3. Lower Bound of the Mean Squared Error of a Biased Estimator

Formally, an estimator is a measurable function that takes as input a random variable (often called observations), and with an output space that contains the vector of parameters to estimate. In the estimation theory, it is known that among all estimators, with the same input random variable and the same output space, there exists a lower bound for the mean squared error [22]. Therefore, it is interesting to construct it and compare it with the mean squared error of the considered estimator.

The misspecified Cramer-Rao Bound (MCRB) is defined as [23]:

$$\text{MCRB} \triangleq (I_n + \Psi_{\theta}) I_{\theta}^{-1} (I_n + \Psi_{\theta})^T + b[\hat{\theta}] b[\hat{\theta}]^T \qquad (17)$$

where $I_n$ represents the identity matrix of size $n$.

The term $\Psi_{\theta}$ represents the Jacobian matrix of the estimator's bias, denoted $b[\hat{\theta}]$, at the vector $\theta$, when it is seen as a function of the true vector of parameters $\theta$. The previous section has shown that, in general, the considered estimator $\hat{\theta}$ is biased. Therefore, one defines this bias function as:

$$\psi : \mathbb{R}^n \to \mathbb{R}^n \\ x \mapsto S(g(x) - h(\theta_0)) + \theta_0 - x. \qquad (18)$$

The Jacobian matrix of $\psi$ at $\theta$ is given by:

$$\Psi_{\theta} = SG_{\theta} - I_n \qquad (19)$$
where $G_θ$ indicates the Jacobian matrix of $g$ at $θ$.

The term $I_θ$ denotes the Fisher Information matrix [23]. Since $y \sim N(g(θ), R)$, the likelihood function is given by:

$$l_y : \mathbb{R}^n \to \mathbb{R} \\ x \mapsto \frac{1}{(2\pi)^{\frac{m}{2}} \sqrt{|R|}} \exp\left(-\frac{1}{2} \|y - g(x)\|_{R}^2\right). \qquad (20)$$

Under some regularity conditions [22], the Fisher Information matrix is defined as:

$$I_{\theta} \triangleq -\mathbb{E}[d_{\theta}^{2}(\log \circ l_{y})] \qquad (21)$$

where $d_θ^2 f$ corresponds to the Hessian matrix of a function $f$ at $θ$. In the considered case, one can show that:

$$I_{\theta} = G_{\theta}^{\top} R^{-1} G_{\theta}. \qquad (22)$$

Finally, assuming that $G_θ$ represents an injective application, then the Fisher Information matrix invert is well defined.

Finally, a theorem states that [23]:

$$\text{MSE}[\hat{\theta}] \geq \text{MCRB} \qquad (23)$$

where the inequality symbol $\geq$ represents the Loewner partial order.

### 4.4. Application to the Problem

Now that a theoretical framework has been written, it can be used to estimate a GNSS receiver displacement and clock bias difference from a misspecified TDCP observation model and quantify the impact of this misspecified model on the estimated displacement. The link between the different quantities introduced in the previous sections is summarized in Table 1. In the next section, the receiver displacement and clock bias difference will be estimated from a set of TDCP observations of $m$ satellites.

<table><caption>**Table 1.** Link between the terms introduced in Sections 3 and 4.</caption><thead><tr><th>Estimation</th><th>TDCP</th></tr></thead><tbody><tr><td>θ as in Equation (9)</td><td>(r<sub>1,2</sub><sup>T</sup>, dt<sub>1,2</sub><sup>r</sup>)<sup>T</sup> as in Equation (4)</td></tr><tr><td>y as in Equation (9)</td><td>λ((Φ<sub>2</sub> - Φ<sub>1</sub>)<sup>1</sup>, ..., (Φ<sub>2</sub> - Φ<sub>1</sub>)<sup>m</sup>)<sup>T</sup><br>with λ(Φ<sub>2</sub> - Φ<sub>1</sub>)<sup>i</sup> as in Equation (3)</td></tr><tr><td>g as in Equation (9)</td><td>(f<sub>r<sub>1</sub></sub><sup>1</sup>, ..., f<sub>r<sub>1</sub></sub><sup>m</sup>), with f<sub>r<sub>1</sub></sub><sup>i</sup> as in Equation (5)</td></tr><tr><td>h as in Equation (10)</td><td>(f<sub>r<sub>1</sub></sub><sup>1</sup>, ..., f<sub>r<sub>1</sub></sub><sup>m</sup>), with f<sub>r<sub>1</sub></sub><sup>i</sup> as in Equation (6)</td></tr><tr><td>R as in Equation (9)</td><td>diag(σ<sub>TDCP<sub>1</sub></sub><sup>2</sup>, ..., σ<sub>TDCP<sub>m</sub></sub><sup>2</sup>)<br>where σ<sub>TDCP<sub>i</sub></sub><sup>2</sup> is the variance of ε Equation (3)</td></tr></tbody></table>

## 5. Experiments

In this experiment, the set of all visible GPS satellites in Toulouse, France, on the first of January 2024, 12:00 UTC are considered. A skyplot is provided in Figure 2. A synthetic GNSS receiver is constructed, with known coordinates $r_1$. Its displacement vector is a sample of a random zero-mean Gaussian random variable with a covariance matrix proportional to the time interval between successive epochs, i.e., $r_{1,2}$ is a sample of
the random variable $N(0_3, (t_2 - t_1)I_3)$, as this stochastic model captures small random movements of the receiver between epochs while maintaining a realistic representation of a pedestrian receiver over short timescales [24,25]. The receiver clock bias difference $dt_{1,2}^r$ is a sample of the law $N(0, t_2 - t_1)$.

The satellite positions are derived from precise orbit data provided in the Standard Product 3 (SP3) format. These SP3 files offer highly accurate satellite ephemerides, providing the real positions of GPS satellites. This ensures that the data reflect actual satellite orbits, further enhancing the realism of the simulation.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Skyplot of GPS satellites in Toulouse, France, on the first of January 2024, 12:00 UTC.</figcaption></figure>

## 5.1. *Experimental Uncertainty Quantification of TDCP Model*

The first experiment characterizes the TDCP model Equation (7), for a single satellite. A satellite is chosen (G17). For each value of $\Delta r_1$, a value for $\Delta f$ Equation (7) is obtained, indicating the discrepancy between the two TDCP models.

In Figures 3 and 4, the set of values taken by $\Delta r_1$ is set to the unit sphere in order to study the importance of its direction. The time between the two epochs is set to $t_2 - t_1 = 1$ s, which is a common value in standard receivers. The figures indicate that the discrepancy between the two models is maximum and minimum on two antipodal points, is null on the equator defined by these two antipodal points, and is continuous, as there is no rapid change of value.

In Appendix A, one shows that the maximum and minimum errors are along the direction given by the vector $s_2 - s_1 - r_{1,2}$, which approximately corresponds to the satellite displacement, and the error is null for orthogonal vectors to this direction.

Overall, the error is low compared to the norm of the first position error, with a maximum error that is approximately $10^4$ lower than $||\Delta r_1||$.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** $||\Delta r_1|| = 1$ m, $t_2 - t_1 = 1$ s. $\mathbf{u}_x$, $\mathbf{u}_y$, and $\mathbf{u}_z$ indicate the x, y, and z components, in the ECEF frame, of $\Delta r_1$. The dot indicates the satellite displacement direction.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** $\|\Delta\mathbf{r}_1\| = 1 \text{ m}, t_2 - t_1 = 1 \text{ s. } \theta \text{ and } \phi \text{ are the spherical coordinates in the ECEF frame of } \Delta\mathbf{r}_1$. The cross indicates the satellite displacement direction.</figcaption></figure>

The same experiment is conducted with a much higher first fix error (i.e., a larger sphere). The results are visible in Figures 5 and 6. Again, the worst direction is given by the vector $\mathbf{s}_2 - \mathbf{s}_1 - \mathbf{r}_{1,2}$, and one observes a ratio between $\|\Delta\mathbf{r}_1\|$ and $\Delta f$ approximately equal to $10^4$ in the worst cases.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** $\|\Delta\mathbf{r}_1\| = 1 \text{ km}, t_2 - t_1 = 1 \text{ s. } \mathbf{u}_x, \mathbf{u}_y$, and $\mathbf{u}_z$ indicate the $x, y$, and $z$ components, in the ECEF frame, of $\frac{\Delta\mathbf{r}_1}{\|\Delta\mathbf{r}_1\|}$. The dot indicates the satellite displacement direction.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** $\|\Delta\mathbf{r}_1\| = 1 \text{ km}, t_2 - t_1 = 1 \text{ s. } \theta \text{ and } \phi \text{ are the spherical coordinates in the ECEF frame of } \Delta\mathbf{r}_1$. The cross indicates the satellite displacement direction.</figcaption></figure>

The study of the worst case direction is further extended, this time by observing the values taken by $\Delta f$ on the set $\Delta\mathbf{r}_1 \propto \mathbf{s}_2 - \mathbf{s}_1 - \mathbf{r}_{1,2}$, for multiple time intervals $t_2 - t_1$. The results are visible in Figure 7. For convenience, a log-log plot is used. For all the considered time intervals $t_2 - t_1$, the plots indicate a linear relation between the model mismatch $\Delta f$ and the norm of the first-fix error $\|\Delta\mathbf{r}_1\|$ since the slope of all the curves is equal to 1 and from Equation (8), $\Delta f(0) = 0$. For each value of the time interval, the linear coefficient can be read at $\|\Delta\mathbf{r}_1\| = 1 \text{ m}$. For example, in the case of $t_2 - t_1 = 1 \text{ s}$, this coefficient is equal to $1.2 \times 10^{-4}$. See Appendix B for a computation of these values. It is visible that the linear coefficient depends on the time interval: the greater the time interval, the greater the
coefficient. One can conclude that for a given value of ||Δr₁||, the greater the time interval, the greater the model mismatch.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Δr₁ ~ s₂ - s₁ - r₁,₂, which corresponds to the worst direction.</figcaption></figure>

Finally, the largest value taken by Δf is computed, for the product of a set of values of ||Δr₁|| and of a set of time intervals between the two observation epochs. This is done in order to explore how the movement of satellites over time would affect the model mismatch. By extending the time between successive observations, the GNSS satellites in the sky change their positions more significantly due to their dynamics. The purpose of this application is to capture a wider range of applications, which can influence the accuracy and reliability of the positioning solutions. These values are visible in Table 2. Omitting rounding errors, the values taken show a bilinear relation with ||Δr₁|| and t₂ - t₁. This result is hard to predict, as Δf is not a linear function of Δr₁, and is even harder to predict for t₂ - t₁, as the time interval does not explicitly appear in Δf.

<table><caption>**Table 2.** Maximum value of Δf Equation (7) with respect to ||Δr₁|| and t₂ - t₁.</caption><thead><tr><th>t<sub>2</sub> - t<sub>1</sub><br>||Δr<sub>1</sub>||</th><th>1 s</th><th>1 min</th><th>30 min</th></tr></thead><tbody><tr><td>1 m</td><td>1.2 × 10<sup>-4</sup> m</td><td>7.1 × 10<sup>-3</sup> m</td><td>2.1 × 10<sup>-1</sup> m</td></tr><tr><td>10 m</td><td>1.2 × 10<sup>-3</sup> m</td><td>7.1 × 10<sup>-2</sup> m</td><td>2.1 m</td></tr><tr><td>1 km</td><td>1.2 × 10<sup>-1</sup> m</td><td>7.1 m</td><td>2.1 × 10<sup>1</sup> m</td></tr><tr><td>1000 km</td><td>1.2 × 10<sup>2</sup> m</td><td>7.1 × 10<sup>3</sup> m</td><td>2.1 × 10<sup>5</sup> m</td></tr></tbody></table>

Now that the model mismatch has been characterized for a single satellite, multiple satellites are going to be used to estimate the receiver's displacement and clock bias difference.

## 5.2. Experimental Uncertainty Quantification of Estimated Displacement with Noiseless Measurements

In this experiment, the receiver's true positions r₁ and r₂ are set as previously explained, and the receiver's clock bias difference follows a zero-mean Gaussian random variable, with a variance proportional to the time interval t₂ - t₁. The set of all visible GPS satellites at the two epochs is used to generate noiseless synthetic TDCP measurements, as in Equation (3), with σ<sub>TDCP</sub> = 0. They are used to estimate the receiver displacement and clock bias difference between the two epochs, following Equation (6). The choice regarding noiseless measurements is made in order to have an error of displacement estimation due to a first-fix error only.

For the metric, since the interest is in the estimation of the receiver displacement vector r̂<sub>1,2</sub>, the mean squared error matrix is expressed in a basis such that the submatrix corresponding to the geometrical part is extracted and is denoted MSE, and the square root of the trace of this matrix is computed, i.e., √Tr(MSE). This choice is made in order
to plot the mean Euclidean distance error. Also, denoting PDOP the position dilution of precision, this metric is equal to $\sigma_{\text{TDCP}} \times \text{PDOP}$, in the case of independent and identically distributed noise [26].

In Figures 8 and 9, the set of values taken by $\Delta\mathbf{r}_1$ corresponds to two spheres with different radii, i.e., for two values of $\|\Delta\mathbf{r}_1\|$. The time interval is set to $t_2 - t_1 = 1$ s. These two plots have the same shape: two modes on the sphere, whose maximum values are close to antipodal and no direction where the estimation error is null. Although these figures look similar, they differ by the magnitude of the values. The local Up and North directions are displayed, but they do not show any link with the errors. More generally, experiments lead to the conclusion that for any value of $\|\Delta\mathbf{r}_1\|$, the obtained plot has a shape similar to the ones in Figures 8 and 9, and that the set of values taken by $\sqrt{\text{Tr}(\text{MSE})}$ is linear with respect to $\|\Delta\mathbf{r}_1\|$.

In order to show these conclusions, the estimation error is studied in a given direction, which is defined by the angles $\theta = \frac{\pi}{2}$ and $\phi = \frac{\pi}{2}$ in the ECEF frame. The reason why this direction is chosen is that it corresponds approximately to the direction where $\sqrt{\text{Tr}(\text{MSE})}$ is maximum for a given value of $\|\Delta\mathbf{r}_1\|$.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** $\sqrt{\text{Tr}(\text{MSE})}$ for noiseless measurements, with $\|\Delta\mathbf{r}_1\| = 1$ m, $t_2 - t_1 = 1$ s. $\theta$ and $\phi$ are the spherical coordinates, in the ECEF frame, of $\Delta\mathbf{r}_1$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** $\sqrt{\text{Tr}(\text{MSE})}$ for noiseless measurements, with $\|\Delta\mathbf{r}_1\| = 1$ km, $t_2 - t_1 = 1$ s. $\theta$ and $\phi$ are the spherical coordinates, in the ECEF frame, of $\Delta\mathbf{r}_1$.</figcaption></figure>

First, in Figure 10, the estimation error is displayed as a function of $\|\Delta\mathbf{r}_1\|$ for multiple time intervals. It shows a linear relation between the estimation error and the norm of the first-fix error for all the time intervals, as the slopes of the curves are all equal to 1. The linear coefficient between $\sqrt{\text{Tr}(\text{MSE})}$ and $\|\Delta\mathbf{r}_1\|$ can be read for $\|\Delta\mathbf{r}_1\| = 1$ m: the greater the time interval, the greater the coefficient. Again, the linearity is hard to predict, as writing the value of $\sqrt{\text{Tr}(\text{MSE})}$ as a linear function of $\Delta\mathbf{r}_1$, when restricted to a specific direction, is not straightforward.

In Figure 11, the estimation error is displayed as a function of $t_2 - t_1$, for multiple values of $\|\Delta\mathbf{r}_1\|$. Again, it shows a linear relation between the estimation error $\sqrt{\text{Tr}(\text{MSE})}$
and the time interval $t_2 - t_1$, for all the values of $||\Delta r_1||$. The linear coefficient between $\sqrt{\text{Tr}(\text{MSE})}$ and $t_2 - t_1$ can be read for $t_2 - t_1 = 1$ s: the greater is $||\Delta r_1||$, the greater is the coefficient. This result is hard to predict, as $\sqrt{\text{Tr}(\text{MSE})}$ is not an explicit function of $t_2 - t_1$.

Finally, Tables 3 and 4 present visible minimum and maximum values of $\sqrt{\text{Tr}(\text{MSE})}$, respectively, for multiple values of $t_2 - t_1$ and for multiple values of $||\Delta r_1||$. Decomposing $\sqrt{\text{Tr}(\text{MSE})}$ in its horizontal and vertical components, in Tables 5 and 6 are the maximum value of $\sqrt{\text{Tr}(\text{MSE})}$ in the horizontal plane and vertical line, respectively, for multiple values of $t_2 - t_1$ and for multiple values of $||\Delta r_1||$. Omitting rounding errors and cases where both $||\Delta r_1||$ and $t_2 - t_1$ are very large, both the minimum and maximum values are linear with respect to $||\Delta r_1||$, but not linear with respect to $t_2 - t_1$. The non-linearity with respect to time is due to the disappearance of G14 after a few minutes, meaning it is not considered in the solution for $t_2 - t_1 = 30$ min. Again, the linearity is hard to predict, as writing the maximum or minimum values of $\sqrt{\text{Tr}(\text{MSE})}$ on a sphere, as a linear function of $||\Delta r_1||$ is not straightforward. From these tables, conclusions can be made regarding the validity of the first-fix approximation. For example, if one has to estimate a receiver displacement such that $||\hat{r}_{1,2} - r_{1,2}|| \le 1.3$ mm, for a time interval of 1 s, for an initial error of the worst amplitude $||\Delta r_1|| = 10$ m, then the approximation $r_1 \approx \hat{r}_1$ is not acceptable, without considering any other source of perturbation.

Overall, it is experimentally shown that the greater $||\Delta r_1||$, the greater $\sqrt{\text{Tr}(\text{MSE})}$, with a linear factor between these two quantities that depends on the time interval $t_2 - t_1$.

<table><caption>**Table 3.** Min values of $\sqrt{\text{Tr}(\text{MSE})}$ with respect to $||\Delta r_1||$ and $t_2 - t_1$.</caption><thead><tr><th>$t_2 - t_1$</th><th>1 s</th><th>1 min</th><th>30 min</th></tr></thead><tbody><tr><th>1 m</th><td>2.4 × 10<sup>-5</sup> m</td><td>1.3 × 10<sup>-3</sup> m</td><td>7.7 × 10<sup>-2</sup> m</td></tr><tr><th>10 m</th><td>2.4 × 10<sup>-4</sup> m</td><td>1.3 × 10<sup>-2</sup> m</td><td>7.7 × 10<sup>-1</sup> m</td></tr><tr><th>1 km</th><td>2.4 × 10<sup>-2</sup> m</td><td>1.3 m</td><td>7.7 × 10<sup>1</sup> m</td></tr><tr><th>1000 km</th><td>2.4 × 10<sup>1</sup> m</td><td>1.3 × 10<sup>3</sup> m</td><td>7 × 10<sup>4</sup> m</td></tr></tbody></table>

<table><caption>**Table 4.** Max values of $\sqrt{\text{Tr}(\text{MSE})}$ with respect to $||\Delta r_1||$ and $t_2 - t_1$.</caption><thead><tr><th>$t_2 - t_1$</th><th>1 s</th><th>1 min</th><th>30 min</th></tr></thead><tbody><tr><th>1 m</th><td>1.3 × 10<sup>-4</sup> m</td><td>7.7 × 10<sup>-3</sup> m</td><td>3.8 × 10<sup>-1</sup> m</td></tr><tr><th>10 m</th><td>1.3 × 10<sup>-3</sup> m</td><td>7.7 × 10<sup>-2</sup> m</td><td>3.8 m</td></tr><tr><th>1 km</th><td>1.3 × 10<sup>-1</sup> m</td><td>7.7 m</td><td>3.8 × 10<sup>2</sup> m</td></tr><tr><th>1000 km</th><td>1.3 × 10<sup>2</sup> m</td><td>7.8 × 10<sup>3</sup> m</td><td>3.8 × 10<sup>5</sup> m</td></tr></tbody></table>

<table><caption>**Table 5.** Max values of $\sqrt{\text{Tr}(\text{MSE})}$ in the receiver local horizontal plane with respect to $||\Delta r_1||$ and $t_2 - t_1$.</caption><thead><tr><th>$t_2 - t_1$</th><th>1 s</th><th>1 min</th><th>30 min</th></tr></thead><tbody><tr><th>1 m</th><td>8.8 × 10<sup>-5</sup> m</td><td>5.3 × 10<sup>-3</sup> m</td><td>1.7 × 10<sup>-1</sup> m</td></tr><tr><th>10 m</th><td>8.8 × 10<sup>-4</sup> m</td><td>5.3 × 10<sup>-2</sup> m</td><td>1.7 m</td></tr><tr><th>1 km</th><td>8.8 × 10<sup>-2</sup> m</td><td>5.3 m</td><td>1.7 × 10<sup>2</sup> m</td></tr><tr><th>1000 km</th><td>8.9 × 10<sup>1</sup> m</td><td>5.3 × 10<sup>3</sup> m</td><td>1.7 × 10<sup>5</sup> m</td></tr></tbody></table>
<table><caption>**Table 6.** Max values of $\sqrt{\text{Tr}(\text{MSE})}$ on the receiver local vertical line with respect to $\|\Delta\mathbf{r}_1\|$ and $t_2 - t_1$.</caption><thead><tr><th>$\|\Delta\mathbf{r}_1\|$</th><th>1 s</th><th>1 min</th><th>30 min</th></tr></thead><tbody><tr><td>1 m</td><td>1.1 × 10<sup>-4</sup> m</td><td>6.7 × 10<sup>-3</sup> m</td><td>3.2 × 10<sup>-1</sup> m</td></tr><tr><td>10 m</td><td>1.1 × 10<sup>-3</sup> m</td><td>6.7 × 10<sup>-2</sup> m</td><td>3.2 m</td></tr><tr><td>1 km</td><td>1.1 × 10<sup>-1</sup> m</td><td>6.7 m</td><td>3.2 × 10<sup>2</sup> m</td></tr><tr><td>1000 km</td><td>1.2 × 10<sup>2</sup> m</td><td>6.7 × 10<sup>3</sup> m</td><td>3.2 × 10<sup>5</sup> m</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** $\sqrt{\text{Tr}(\text{MSE})}$ for noiseless measurements and for multiple values of $t_2 - t_1$ as a function of $\|\Delta\mathbf{r}_1\|$, with $\Delta\mathbf{r}_1$ in the direction given by the angles $\theta = \frac{\pi}{2}$ and $\phi = \frac{\pi}{2}$ in the ECEF frame.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** $\sqrt{\text{Tr}(\text{MSE})}$ for noiseless measurements and for multiple values of $\|\Delta\mathbf{r}_1\|$, as a function of $t_2 - t_1$, with $\Delta\mathbf{r}_1$ in the direction given by the angles $\theta = \frac{\pi}{2}$ and $\phi = \frac{\pi}{2}$ in the ECEF frame.</figcaption></figure>

### 5.3. Experimental Uncertainty Quantification of Estimated Displacement with Noisy Measurements

In the same scenario as the previous section, noisy measurements are generated following Equation (3), with $\sigma_{\text{TDCP}} = \sqrt{2}\sigma_{\phi}$ and $\sigma_{\phi} = \lambda_{L1} \times 10^{-1}$ m, with $\lambda_{L1}$ denoting the wavelength of the L1 frequency band. The time interval is set to $t_2 - t_1 = 1$ s. Note that alternative noise models could also be considered in order to reflect more realistic weighting schemes. For instance, by incorporating satellite elevation-dependent variances.

As explained previously, since the interest is in the estimation of the receiver displacement vector $\hat{\mathbf{r}}_{1,2}$, the quantity $\sqrt{\text{Tr}(\text{MSE})}$ is computed by extracting the geometrical part of the MSE matrix. The same process is applied to the MCRB matrix, and the square root of the trace is computed. This choice is justified by the fact that the Loewner partial order has the following property: for all pairs of real symmetric positive semi-definite matrices $A$ and $B$, $A \ge B \Rightarrow \sqrt{\text{Tr}(A)} \ge \sqrt{\text{Tr}(B)}$ (in particular, $\sqrt{\text{Tr}(A)} \ge 0$).

Figure 12 shows $\sqrt{\text{Tr}(\text{MSE})}$ for $\|\Delta\mathbf{r}_1\| = 1$ m. On the same set, the difference between $\sqrt{\text{Tr}(\text{MSE})}$ and its lower bound $\sqrt{\text{Tr}(\text{MCRB})}$ is visible in Figure 13. The results indicate that, as was the case in the absence of noise, the error is characterized by two modes on
the sphere. Also, the estimator's error is close to the lowest value that it can reach, with no particular direction appearing.

Figures 14 and 15 depict the same experiment, with $||\Delta\mathbf{r}_1|| = 1$ km. Again, the results indicate that the estimator's error is close to the lowest value that it can reach.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** $\sqrt{\text{Tr}(\text{MSE})}$ for noisy measurements, with $||\Delta\mathbf{r}_1|| = 1$ m, $t_2 - t_1 = 1$ s. $\theta$ and $\phi$ are the spherical coordinates in the ECEF frame of $\Delta\mathbf{r}_1$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Difference between $\sqrt{\text{Tr}(\text{MSE})}$ and $\sqrt{\text{Tr}(\text{MCRB})}$ for noisy measurements with $||\Delta\mathbf{r}_1|| = 1$ m, $t_2 - t_1 = 1$ s. $\theta$ and $\phi$ are the spherical coordinates in the ECEF frame of $\Delta\mathbf{r}_1$. The negative values exist due to numerical errors.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** $\sqrt{\text{Tr}(\text{MSE})}$ for noisy measurements, with $||\Delta\mathbf{r}_1|| = 1$ km, $t_2 - t_1 = 1$ s. $\theta$ and $\phi$ are the spherical coordinates in the ECEF frame of $\Delta\mathbf{r}_1$.</figcaption></figure>

Finally, as done previously, the estimation error is studied along an axis defined by the angles $\theta = \frac{\pi}{2}$ and $\phi = \frac{\pi}{2}$ in the ECEF frame. Figure 16 shows the error and its lower bound. At $\Delta\mathbf{r}_1 = 0$, the estimation error is due to the variance of the additive noise within the observations. This noise term adds up to the bias as $||\Delta\mathbf{r}_1||$ increases, which produces the curve that is observed. For large values of $||\Delta\mathbf{r}_1||$, the error is mainly due to the misspecified model. Also, the results indicate that the estimator's error has reached the lowest value possible.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Difference between $\sqrt{\text{Tr}(\text{MSE})}$ and $\sqrt{\text{Tr}(\text{MCRB})}$ for noisy measurements with $\|\Delta\mathbf{r}_1\| = 1$ km, $t_2 - t_1 = 1$ s. $\theta$ and $\phi$ are the spherical coordinates in the ECEF frame of $\Delta\mathbf{r}_1$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** $\sqrt{\text{Tr}(\text{MSE})}$ and $\sqrt{\text{Tr}(\text{MCRB})}$ for noisy measurements as a function of $\|\Delta\mathbf{r}_1\|$, with $t_2 - t_1 = 1$ s and $\Delta\mathbf{r}_1$ in the direction given by the angles $\theta = \frac{\pi}{2}$ and $\phi = \frac{\pi}{2}$ in the ECEF frame.</figcaption></figure>

In the last experiment, Figure 17 depicts the maximum value of $\sqrt{\text{Tr}(\text{MSE})}$ as a function of $\|\Delta\mathbf{r}_1\|$.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Maximum value of $\sqrt{\text{Tr}(\text{MSE})}$ and corresponding $\sqrt{\text{Tr}(\text{MCRB})}$ for noisy measurements as a function of $\|\Delta\mathbf{r}_1\|$ with $t_2 - t_1 = 1$ s.</figcaption></figure>

To conclude this section, one can say that all the estimators built in these experiments are close to the MCRB. The mean error depends on the amount of noise in the observations, and on the misspecification of the model used.

# 6. Conclusions and Future Work

The TDCP model is expressed as a function of a GNSS receiver displacement and clock bias difference. The impact on the estimated displacement of the uncertainty regarding the receiver coordinates at the first epoch has been quantified for a given geometry at a given location. A similar study can be run in real-time by any user anywhere in the world, as the proposed method does not rely on final products, and real-time ephemerides would provide similar results. Considering the displacement estimation error due to a misspecified
first fix only from a set of TDCP observations and given criteria on the precision of the estimated displacement, the provided tables are an indicator of the validity of the first fix. It is experimentally shown that for long time intervals, the impact is important. This is due to the fact that for long time intervals, the satellite coordinates change more, and lead to different receiver-to-satellite unit vectors. This effect would be more prominent when the considered satellites are at a low altitude. Therefore, when considering Low Earth Orbit satellites, whose displacements are much faster than GNSS satellites in Medium Earth Orbits, the impact would be greater. One might also notice that the same experiment could be transposed to an RTK framework: at a given epoch, the station coordinates are known, the baseline vector corresponds to the receiver displacement, and the single difference between receivers would correspond to TDCPs.

**Author Contributions:** Conceptualization: all authors; methodology: all authors; software: H.C.; validation: J.L., J.S. and P.T.; formal analysis: all authors; investigation: all authors; resources: all authors; data curation: all authors; writing—original draft preparation: H.C.; writing—review and editing: J.L., J.S. and P.T.; visualization: all authors; supervision: all authors; project administration: all authors; funding acquisition: J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the EBCON Project under Grant PID2020119244GB-100, and in part by the Consolidated Research Group RAIG under Grant 2021 SGR 00510 of the Departament de Recerca i Universitats de la Generalitat de Catalunya.

**Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

**Conflicts of Interest:** The authors declare that they have no conflicts of interest.

## Appendix A. Direction of Null and Extremum Model Error

Recalling that $\Delta r_1 = \hat{r}_1 - r_1$ and that Equation (8) is the function $\Delta f$ that maps an error $\Delta r_1$ to a model error:

$$ \Delta f(\Delta r_1) = \|s_2 - r_1 - r_{1,2} - \Delta r_1\| \\ - \|s_2 - r_1 - r_{1,2}\| \\ - (\|s_1 - r_1 - \Delta r_1\| - \|s_1 - r_1\|) $$

where $\|x\| = \sqrt{\langle x, x \rangle}$ and $\langle x, y \rangle = x^T y$ is the dot product on $\mathbb{R}^3$. Defining $a \triangleq s_2 - r_1 - r_{1,2}$ and $b \triangleq s_1 - r_1$, one obtains

$$ \Delta f(\Delta r_1) = \|a - \Delta r_1\| - \|a\| - (\|b - \Delta r_1\| - \|b\|). $$

Given that $\Delta f(0) = 0$, one obtains

$$ \Delta f(\Delta r_1) = d(\Delta f)_0(\Delta r_1) + o(\|\Delta r_1\|) \\ = -\left\langle \frac{a}{\|a\|} - \frac{b}{\|b\|}, \Delta r_1 \right\rangle + o(\|\Delta r_1\|) $$

where $d(\Delta f)_0(\Delta r_1)$ is the differential of $\Delta f$ at 0, evaluated at $\Delta r_1$. If this term had to be expressed in the canonical bases of $\mathbb{R}^3$ and $\mathbb{R}$, it would be equal to the Jacobian matrix of $\Delta f$ at 0 multiplied by $\Delta r_1$. Therefore, in a neighborhood of 0, the linear approximation of $\Delta f$ is given by:

$$ \Delta f(\Delta r_1) \approx -\left\langle \frac{a}{\|a\|} - \frac{b}{\|b\|}, \Delta r_1 \right\rangle. $$
Now, one has to study the vector $\frac{a}{\|a\|} - \frac{b}{\|b\|}$ and give it a geometrical interpretation, other than being the difference between two very close unit vectors. One defines the function:

$$ \gamma : \mathbb{R}^3 \to \mathbb{R}^3, x \mapsto \frac{a}{\|a\|} - \frac{a+x}{\|a+x\|}, $$

and defining $v \triangleq b - a = -(s_2 - s_1 - r_{1,2})$, one obtains

$$ \frac{a}{\|a\|} - \frac{b}{\|b\|} = \frac{a}{\|a\|} - \frac{a+v}{\|a+v\|} = \gamma(v). $$

Now that the vector $\frac{a}{\|a\|} - \frac{b}{\|b\|}$ is expressed as a valuation of the function $\gamma$, a linear approximation of it is given: it can be shown that $\forall \alpha \in \mathbb{R}^3 \setminus \{-a\}, \forall x \in \mathbb{R}^3$:

$$ \gamma(\alpha + x) = \gamma(\alpha) + d\gamma_{\alpha}(x) + o(\|x\|) $$

where $d\gamma_{\alpha}(x)$ is the differential of $\gamma$ at $\alpha$, evaluated at $x$. If this term had to be expressed in the canonical basis of $\mathbb{R}^3$, it would be equal to the Jacobian matrix of $\gamma$ at $\alpha$ multiplied by $x$. One can show that:

$$ d\gamma_{\alpha}(x) = - \frac{1}{\|a + \alpha\|} \left( x - \frac{a + \alpha}{\|a + \alpha\|} \left\langle \frac{a + \alpha}{\|a + \alpha\|}, x \right\rangle \right). $$

Since $\gamma(0) = 0$, the linear approximation of the function $\gamma$ around 0 is given by:

$$ \gamma(x) \approx - \frac{1}{\|a\|} \left( x - \frac{a}{\|a\|} \left\langle \frac{a}{\|a\|}, x \right\rangle \right). $$

As the vector $v$ is an infinitesimal and orthogonal vector, with respect to the vector $a = s_2 - r_2$, the scalar product term $\langle \frac{a}{\|a\|}, v \rangle$ is neglected:

$$ \Rightarrow \frac{a}{\|a\|} - \frac{b}{\|b\|} = \gamma(v) \approx - \frac{v}{\|a\|} = \frac{a-b}{\|a\|}. $$

Finally, $\Delta f(\Delta r_1) \approx 0 \Leftrightarrow \langle \frac{a-b}{\|a\|}, \Delta r_1 \rangle \approx 0$, i.e., the model error is null for any error $\Delta r_1$ orthogonal to $a-b = s_2 - s_1 - r_{1,2}$, and the directions of extremum error are given by the vector $s_2 - s_1 - r_{1,2}$.

## Appendix B. Model Error Along the Worst Direction

As a reminder, in Appendix A, the vectors $a \triangleq s_2 - r_1 - r_{1,2}$ and $b \triangleq s_1 - r_1$ are defined, and the vector $v = a - b$ is shown to be the vector that defines the direction of extremum error. Here, the model error $\Delta f$ of Equation (7) is studied along a one-dimensional subspace: for any unit vector $e$, considering the application:

$$ \phi : \mathbb{R} \to \mathbb{R} $$

$$ \lambda \mapsto \Delta f(\lambda e) $$

with

$$ \Delta f(\lambda e) = \|a - \lambda e\| - \|a\| - (\|b - \lambda e\| - \|b\|), $$
one finds out that the derivative of $\phi$ at $\lambda$ is given by:

$$ \phi'(\lambda) = - \frac{\langle \mathbf{a} - \lambda \mathbf{e}, \mathbf{e} \rangle}{\| \mathbf{a} - \lambda \mathbf{e} \|} + \frac{\langle \mathbf{b} - \lambda \mathbf{e}, \mathbf{e} \rangle}{\| \mathbf{b} - \lambda \mathbf{e} \|} $$

Finally, if $e = \frac{\mathbf{v}}{\|\mathbf{v}\|}$ and $\Delta\mathbf{r}_1 = \lambda\mathbf{e}, \lambda \in \mathbb{R}$, then $\|\Delta\mathbf{r}_1\| = |\lambda|$ and one has obtained the derivative of $\phi$ at $\|\Delta\mathbf{r}_1\|$ along the direction of $\mathbf{v}$.

## References

1. Rustamov, R.B.; Hashimov, A.M. Multifunctional Operation and Application of GPS; IntechOpen: Rijeka, Croatia, 2018.

2. Congram, B.; Barfoot, T.D. Experimental Comparison of Visual and Single-Receiver GPS Odometry. arXiv **2021**, arXiv:2106.02122.

3. Kim, J.; Park, M.; Bae, Y.; Kim, O.J.; Kim, D.; Kim, B.; Kee, C. A Low-Cost, High-Precision Vehicle Navigation System for Deep Urban Multipath Environment Using TDCP Measurements. Sensors **2020**, 20, 3254. [CrossRef] [PubMed]

4. Chen, K.; Chang, G.; Chen, C.; Zhu, T. An improved TDCP-GNSS/INS integration scheme considering small cycle slip for low-cost land vehicular applications. Meas. Sci. Technol. **2021**, 32, 055006. [CrossRef]

5. Ji, L.; Sun, R.; Cheng, Q.; Wang, J. Evaluation of the performance of GNSS-based velocity estimation algorithms. Satell. Navig. **2022**, 3, 18. [CrossRef]

6. Sun, R.; Cheng, Q.; Wang, J. Precise vehicle dynamic heading and pitch angle estimation using time-differenced measurements from a single GNSS antenna. GPS Solut. **2020**, 24, 84. [CrossRef]

7. Li, Z.; Wang, D.; Dong, Y.; Wu, J. An enhanced tightly-coupled integrated navigation approach using phase-derived position increment (PDPI) measurement. Optik **2018**, 156, 135–147. [CrossRef]

8. Suzuki, T. GNSS Odometry: Precise Trajectory Estimation Based on Carrier Phase Cycle Slip Estimation. IEEE Robot. Autom. Lett. **2022**, 7, 7319–7326. [CrossRef]

9. Freda, P.; Angrisano, A.; Gaglione, S.; Troisi, S. Time-differenced carrier phases technique for precise GNSS velocity estimation. GPS Solut. **2015**, 19, 335–341. [CrossRef]

10. Lyu, P.; Bai, S.; Lai, J.; Wang, B.; Sun, X.; Huang, K. Optimal Time Difference-Based TDCP-GPS/IMU Navigation Using Graph Optimization. IEEE Trans. Instrum. Meas. **2021**, 70, 1–10. [CrossRef]

11. Olivier, A.; Giovanis, D.G.; Aakash, B.; Chauhan, M.; Vandanapu, L.; Shields, M.D. UQpy: A general purpose Python package and development environment for uncertainty quantification. J. Comput. Sci. **2020**, 47, 101204. [CrossRef]

12. Purnell, D.; Gomez, N.; Chan, N.H.; Strandberg, J.; Holland, D.M.; Hobiger, T. Quantifying the Uncertainty in Ground-Based GNSS-Reflectometry Sea Level Measurements. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. **2020**, 13, 4419–4428. [CrossRef]

13. Martin, A.; Soundy, A.; Panckhurst, B.; Brown, C.; Schumayer, D.; Molteno, T.; Parry, M. Real-time uncertainty quantification using correlated noise models for GNSS positioning. In Proceedings of the 2017 IEEE SENSORS, Glasgow, UK, 29 October–1 November 2017; pp. 1–3.

14. Su, J.; Schön, S. Deterministic Approaches for Bounding GNSS Uncertainty: A Comparative Analysis. In Proceedings of the 10th Workshop on Satellite Navigation Technology (NAVITEC), Noordwijk, The Netherlands, 5–7 April 2022; pp. 1–8.

15. Richmond, C.D.; Horowitz, L.L. Parameter Bounds on Estimation Accuracy Under Model Misspecification. IEEE Trans. Signal Process. **2015**, 63, 2263–2278. [CrossRef]

16. Mennad, A.; Fortunati, S.; El Korso, M.N.; Younsi, A.; Zoubir, A.M.; Renaux, A. Slepian-Bangs-type formulas and the related Misspecified Cramér-Rao Bounds for Complex Elliptically Symmetric distributions. Signal Process. **2018**, 142, 320–329. [CrossRef]

17. Fortunati, S.; Ortega, L. On the efficiency of misspecified Gaussian inference in nonlinear regression: Application to time-delay and Doppler estimation. Signal Process. **2024**, 225, 109614. [CrossRef]

18. Yang, C.; Guo, J.; Geng, T.; Zhao, Q.; Jiang, K.; Xie, X.; Lv, Y. Assessment and Comparison of Broadcast Ionospheric Models: NTCM-BC, BDGIM, and Klobuchar. Remote Sens. **2020**, 12, 1215. [CrossRef]

19. Schüler, T. The TropGrid2 standard tropospheric correction model. GPS Solut. **2014**, 18, 123–131. [CrossRef]

20. Stienne, G.; Reboul, S.; Choquel, J.; Benjelloun, M. Cycle slip detection and repair with a circular on-line change-point detector. Signal Process. **2014**, 100, 51–63. [CrossRef]

21. Nocedal, J.; Wright, S.J. Numerical Optimization; Springer Series in Operations Research and Financial Engineering; Springer: New York, NY, USA, 2006.

22. Nielsen, F. Connected at Infinity II: A Selection of Mathematics by Indians; Chapter Cramér-Rao Lower Bound and Information Geometry; Hindustan Book Agency: Delhi, India, 2013; pp. 18–37.

23. Eldar, Y. Minimum Variance in Biased Estimation: Bounds and Asymptotically Optimal Estimators. Trans. Sig. Proc. **2004**, 52, 1915–1930. [CrossRef]
24. Schöller, C.; Aravantinos, V.; Lay, F.; Knoll, A. What the Constant Velocity Model Can Teach Us About Pedestrian Motion Prediction. *IEEE Robot. Autom. Lett.* **2020**, *5*, 1696–1703. [CrossRef]

25. Pellegrini, S.; Ess, A.; Van Gool, L. *Visual Analysis of Humans: Looking at People*; Chapter Predicting Pedestrian Trajectories; Springer: London, UK, 2011; pp. 473–491.

26. Won, D.H.; Ahn, J.; Lee, S.W.; Lee, J.; Sung, S.; Park, H.W.; Park, J.P.; Lee, Y.J. Weighted DOP with Consideration on Elevation-Dependent Range Errors of GNSS Satellites. *IEEE Trans. Instrum. Meas.* **2012**, *61*, 3241–3250. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.