2.1 Principle of RCM

The process wherein the internal temperature compensated crystal oscillator (TCXO) of the GNSS receiver is replaced with an external oscillator possessing high frequency stability allows modelling its operation in a physically meaningful way over a certain time duration. This method is referred as RCM and it further strengthens the observation geometry over the modelling time duration. In RCM, it is mandatory that the oscillator noise is less than the GNSS receiver noise. Thus, the random integrated frequency fluctuations of the oscillator cannot be resolved by the GNSS observations (code, phase or Doppler).

Oscillator characterisation is typically done based on time-dependent Allan variance $({{\sigma_y}^2} )$ parameter (Allan, Reference Allan1987). This measure allows determining the time prediction errors, which in turn enables to find the time modelling duration over which RCM holds true, i.e. oscillator noise is less than the GNSS receiver noise. Furthermore, the RCM time duration can also be computed based on the power spectral density of the signal (Allan et al., Reference Allan, Hellwig, Kartaschoff, Vanier, Vig, Winkler and Yannoni1988). The noise of the GNSS receiver can be assumed to be about 1% of the signal wavelength in use, i.e. 3 m for L1 C/A code, 0⋅002 m for L1 carrier phase observation. For L1 Doppler observation, the noise can be assumed to be about 0⋅05 m/s.

In the kinematic experiment, a Microsemi MAC SA.35 m oscillator is used. It is a rubidium (Rb) standard atomic clock and its frequency signal serves as clock input to a geodetic grade JAVAD Delta GNSS receiver and a HS low-cost u-blox receiver. To validate the frequency stability information provided in the manufacturer's datasheet and compute the appropriate RCM time duration, a clock characterisation study is carried out in a laboratory at PTB – The National Metrology Institute of Germany, in terms of the short- and long-term frequency stabilities. From this study, it is observed that the estimated frequency stability of MAC SA.35 m agrees with the manufacturer's values. Moreover, the performance of the oscillator is also evaluated in a highly dynamic environment and all the results are shown in Jain et al. (Reference Jain, Krawinkel, Schön and Bauch2020).

The Allan deviation (ADEV) values obtained from the MAC SA.35 m manufacturer's datasheet and computed from laboratory dataset for averaging time from one second to about 17 min are shown Figure 1. In addition, ADEV's for GPS L1 C/A code, phase and Doppler observations for the same averaging time are depicted in Figure 1. As code and carrier-phase observations are phase-based measurements, the noise of C/A code and phase observations are modelled as white phase modulation (WPM) over time, whereas the noise of Doppler is modelled as white frequency modulation (WFM) over time. As shown in Figure 1, the noise of MAC SA.35 m is smaller than the noise of GPS L1 C/A code and Doppler observations but higher than the carrier phase observations. Hence, RCM cannot be applied with phase observations using this oscillator. The estimated ADEV's from laboratory dataset for MAC SA.35 m are transformed to spectral density ${h_\alpha }$ coefficient (Barnes et al., Reference Barnes, Chi, Cutler, Healey, Leeson, McGunigal, Mullen, Smith, Sydnor, Vessot and Winkler1971). Using this ${h_\alpha }$ coefficients specified in Jain et al. (Reference Jain, Krawinkel, Schön and Bauch2020), the spectral behaviour of the oscillator can be modelled through the process noise matrix when estimating state parameters (e.g. position, velocity and timing) through a Kalman filtering approach.

Figure 1. ADEV of microsemi MAC SA.35 m from manufacturer's datasheet and derived from characterization carried out in a laboratory at PTB. GPS L1 C/A code and L1 carrier phase observation noise modelled as white phase modulation (WPM); GPS L1 Doppler noise modelled as white frequency modulation (WFM)

2.2 Kalman filter estimation

In the case of GNSS-based PVT solution, different state errors are estimated largely using a LKF. Moreover, this technique is used in variety of other applications due to its performance, computational efficiency and applicability in real time. It is the best possible linear-estimator with minimum variance when optimum values for measurement and process noise are known (Brown and Hwang, Reference Brown and Hwang2012). In comparison with state estimation with least-square adjustment (LSA), Kalman filter integrates system dynamics along with appropriate process noise models. Different motion models could be used to account for dynamics. During phases when the motion model becomes invalid, the filter can be tuned with process noise up to a certain limit, i.e. uncertainty in the states can be bounded.

The observation equation for pseudo-range measurements in GNSS-based navigation solution between a receiver A and satellite i is

(2.1)\begin{equation}\bar{l}_k^i = \; \sqrt {{{({{X^i} - {X_A}} )}^2} + {{({{Y^i} - {Y_A}} )}^2} + {{({{Z^i} - {Z_A}} )}^2}} + T_A^i\; + I_A^i + c({\delta {t_A} - \delta {t^i}} )+ \delta t_{A,\; rel}^i + \; \epsilon _A^i\end{equation}

where $X,\; Y$ and Z = respective geocentric coordinates for the receiver A and satellite i; $T_A^i$ and $I_A^i$ = tropospheric and ionospheric signal delays; c = speed of light; $\delta {t_A}$ and $\delta {t^i}$ = respective receiver and satellite clock errors; $\delta t_{A,\; rel}^i$ = relativistic time offset and $\epsilon _A^i$ = measurement noise. The ionospheric and tropospheric (wet and dry components) delays are accounted using IGS ionospheric exchange total electron content (IONEX TEC) map (Schaer et al., Reference Schaer, Gurtner and Feltens1998) and the Vienna mapping function 3 (Landskron and Böhm, Reference Landskron and Böhm2018), respectively. In addition, the satellite clock errors and relativistic time offsets are corrected by the orbit and clock products from the multi-GNSS experiment (MGEX) of IGS (Montenbruck et al., Reference Montenbruck, Steigenberger, Prange, Deng, Zhao, Perosanz, Romero, Noll, Stürze, Weber, Schmid, MacLeod and Schaer2017). $\bar{l}_k^i$ is computed for every visible satellite along the complete time series of the dataset. Correspondingly, by using the observation equation for Doppler measurements and appropriate modelling, Doppler observations are computed. Finally, computed code and Doppler observations form the observation vector ${\bar{\textbf{l}}_k}$.

Due to the nonlinear relation between the recorded observations ${\textbf{l}_k}$ and state vector ${\textbf{x}_k}$, the observed-minus-computed (OMC) vector $\Delta {\textbf{l}_k}$ is formed, i.e. the difference between the recorded code and Doppler observations ${\textbf{l}_k}$ and computed observations ${\bar{\textbf{l}}_k}$. With a LKF, the error state vector $\Delta {\hat{\textbf{x}}_k}$ is computed and then added to the a priori state ${\bar{\textbf{x}}_k}$ to obtain the final estimate state ${\hat{\textbf{x}}_k}.$

For the initialisation epoch, the error state and its corresponding variance–covariance matrix (VCM) are set as

(2.2)\begin{equation}\Delta \hat{\textbf{x}}_0^-= 0;\textbf{Q}_{\hat{x}\hat{x},0}^-= {\textbf{Q}_{xx,0}}\end{equation}

Apart from the initialisation step, the error state vector and its VCM are filtered in the update step:

(2.3)\begin{equation}\Delta \hat{\textbf{x}}_k^ += \Delta \hat{\textbf{x}}_k^ -{+} \; {\textbf{K}_k}({\Delta {\textbf{l}_k} - {\textbf{A}_k}\Delta \hat{\textbf{x}}_k^ - } )\; \; \; \textbf{Q}_{xx,k}^ += ({\textbf{I} - {\textbf{K}_k}{\textbf{A}_k}} )\textbf{Q}_{xx,k}^ - {({\textbf{I} - {\textbf{K}_k}{\textbf{A}_k}} )^\textrm{T}} + \; {\textbf{K}_k}{\textbf{Q}_{ll,k}}\textbf{K}_k^\textrm{T}\end{equation}

with the design matrix ${\textbf{A}_k}$ containing the partial derivatives of the functional model evaluated at the nominal trajectory. The observation uncertainty ${\textbf{Q}_{ll,k}}$ is derived using SIGMA- variance model. It was first proposed by Hartinger and Brunner (Reference Hartinger and Brunner1999) for GPS carrier phase observations and is given as

(2.4)\begin{equation}\sigma _{i,\varepsilon }^2 = {C_v}\; {10^{ - \frac{{C/{N_0}}}{{10}}}},\; {\boldsymbol{Q}_{ll}} = \textrm{diag}({\sigma_{i,\varepsilon }^2} )\end{equation}

where $C/{N_0}$ = measured carrier to noise density ratio of the signal in dB-Hz, and ${C_v}$ = variance factor dependent on the receiver–antenna configuration and is estimated from the transformed samples. The applicability of the SIGMA- variance model was further extended and proven with GPS code observations; an improvement in position accuracy was shown to be about 30–50% using real urban data recorded with different types of HS receivers (Wieser et al., Reference Wieser, Gaggl and Hartinger2005). Overall, the model is based on variance factors that are pre-derived for each specific receiver–antenna combination and different GNSS systems based on recorded $C/{N_0}$ of corresponding signals in an open-sky static environment. It must be noted that these calibrated variance factors are not specific to any surrounding in regards to its utilisation.

To compute the variance factors for both code and Doppler observations, a static study is carried out with exactly the same measurement setup. For multiple receivers of the same type, the variance factors obtained after post-processing the static dataset are almost similar considering the respective GNSS system. Hence, its mean value is computed and listed in Table 1 for the code observations. For the Doppler observation variance factors, the reader is referred to Jain et al. (Reference Jain, Kulemann and Schön2021). These variance factors are then used to model the variance of the corresponding pseudo-range and Doppler observations in relation to the measured $C/{N_0}$ during the kinematic experiment.

Table 1. Estimated code observation variance factors for different receiver-antenna combinations

The Kalman gain ${\textbf{K}_k}$ is evaluated as

(2.5)\begin{equation}{\textbf{K}_k} = \textbf{Q}_{xx,k}^ - \textbf{A}_k^\textrm{T}{({{\textbf{Q}_{ll,k}} + {\textbf{A}_k}\textbf{Q}_{xx,k}^ - \textbf{A}_k^\textrm{T}} )^{ - 1}}\end{equation}

Later, the prediction step is performed where the updated error states are time propagated:

(2.6)\begin{equation}\Delta \textbf{x}_{k + 1}^ -= {\mathbf{\Phi }_k}\Delta \textbf{x}_k^ + ;\; \; \textbf{Q}_{xx,k + 1}^ -= {\mathbf{\Phi }_k}\boldsymbol{Q}_{xx,k}^ - \mathbf{\Phi }_k^T + {\textbf{Q}_{\omega\omega,k}}\end{equation}

where ${\mathbf{\Phi }_{\boldsymbol{k}}}$ = transition matrix and ${\textbf{Q}_{\omega\omega,k}}$ = process noise matrix. For the kinematic dataset processing, the transition matrix is selected based on a constant velocity (CV) model. The process noise for the corresponding position and velocity states are modelled as a random ramp and a random walk process, respectively. In the case of RCM, the elements of clock process noise matrix (${\textbf{Q}_{ \omega\omega- \textrm{clk},k}}$) are computed using the spectral density ${h_\alpha }$ coefficients of the oscillator in use and is given as

(2.7)\begin{equation}{\textbf{Q}_{\omega\omega- \textrm{clk},k}} = \left[ {\begin{array}{*{20}{l}} {{h_0}\dfrac{{\Delta t}}{2} + {h_{ - 1}}2\Delta {t^2} + {h_{ - 2}}\dfrac{{2{\pi^2}\Delta {t^3}}}{3}}& {{h_{ - 1}}\Delta t + {h_{ - 2}}{\pi^2}\Delta {t^2}}\\ {{h_{ - 1}}\Delta t + {h_{ - 2}}{\pi^2}\Delta {t^2}}& {{h_0}\dfrac{1}{{2\Delta t}} + 4{h_{ - 1}} + {h_{ - 2}}\dfrac{{8{\pi^2}\Delta t}}{3}} \end{array}} \right]\end{equation}

where ${h_0} = 8 \times {10^{ - 22}},\; \; {h_{ - 1}} = 7 \cdot 2 \times {10^{ - 25}}$ and ${h_{ - 2}} = 3 \times {10^{ - 29}}$ are the spectral density coefficients of white frequency noise (WFN), flicker frequency noise (FFN) and random-walk frequency noise (RWFN) applied in this study, respectively. More information about the computation of the process noise matrix in the absence of RWFN spectral coefficient and its implementation in KF is explained by Krawinkel (Reference Krawinkel2018). For the case when RCM is not applied, the clock time offset is modelled as a random ramp process, whereas the clock frequency offset is modelled as a random walk process. It must be noted that different results of the KF with and without applying RCM are exclusively due to the differences in the clock process noise models.

3.1 Experiment measurement setup

The measurement configuration and the vehicle used during the experiment is shown in Figures 2(a) and 2(b), respectively. It consists of two geodetic-grade JAVAD GNSS receivers [type: Delta TRE-G3 T(H)], and two HS u-blox receivers (type: NEO M8 T). Both JAVAD GNSS receivers were running on the same firmware version, and the same applies for the u-blox receivers. In addition to the four receivers, there was a Microsemi miniaturised atomic clock (MAC) SA.35 m, a Novatel SPAN receiver connected with iMAR inertial measurement unit (IMU) and an active GNSS splitter, which was connected to the geodetic antenna (type: JAVAD GrAnt G3 T) placed at the top of the vehicle.

Figure 2. Kinematic experiment in hannover. (a) measurement setup, (b) institute vehicle used during the experiment

The internal circuitry of JAVAD receiver (0082) and u-blox receiver (1771) were driven by MAC SA.35 m, whereas the other JAVAD receiver (0081) and u-blox receiver (0867) were not modified and were driven by their own internal TCXO. For the u-blox receiver (1771), there is no external frequency input as for the JAVAD receivers, which makes connecting an external clock straightforward. The TCXO within u-blox (1771) was removed; the 10 MHz signal from the MAC SA.35 m was modulated with a signal wave generator (SWG) to 26 MHz (i.e. u-blox clock signal frequency) and fed via a connector to the receiver. The feasibility and PVT performance of this modified u-blox receiver (1771) in an open-sky static environment is explained by Bochkati and Schön (Reference Bochkati and Schön2018). The Novatel SPAN receiver in conjunction with the iMAR IMU forms a multisensor system and records IMU data in GPS time frame. From this IMU data, GNSS data from a reference station located at our institute laboratory terrace and GNSS data from JAVAD receiver (0082), a tightly coupled (TC) reference trajectory solution was computed for the complete duration of the urban experiment. It must be noted that the computation is done in relative carrier phase-based GNSS post-processing mode using the TerraPos software (Kjørsvik et al., Reference Kjørsvik, Øvstedal, Gjevestad and Sideris2009).

3.2 Test drive details

The kinematic experiment was carried out on 22 May 2019 around the western and southern part of the city of Hannover, Germany. The experiment lasted for a duration of about two hours and the complete test route is shown in Figure 3. The test route is comprised of a variety of light and moderate urban conditions, ranging from open-sky to overhanging trees and branches, to narrow roads surrounded by buildings of heights in the range of about 20–30 meters on both sides. Such conditions lead to multipath, blockage or attenuation of signals and interference with other high-intensity signals from the adjacent GNSS frequency bands. The measurement starts with a static phase of about seven minutes in a parking area; later, the vehicle starts moving. It first passes through the inner-city area, which has a more open-sky characteristic. Following this, the vehicle moves along Marienstraße, where the street is mostly flanked on both sides by buildings (Figure 4c). Later, a certain loop, shown in right part of Figure 3, was driven six times, with the first loop being slightly larger and its beginning is marked with point B. The following five loops were exactly similar, taking approximately seven minutes per loop and its ending is marked with point C. Some pictures of the streets along the loop route are shown in Figures 4(a) and 4(b). It consists of mostly three-storey buildings on both sides with some tress and very narrow streets, with the width being in the range of about 4-6 m. Later, another test loop is driven more than two times and its start is marked with point D in Figure 3. This loop route includes Marienstraße, two underpasses and an overhead bridge. Figure 4(d) is an example street view picture of one such underpass. Moreover, there are several standstill phases at traffic signals and crossings along the complete test route. The details of the marked points on the urban trajectory map and their approximate GPS times are listed in Table 2.

Figure 3. Full measurement drive route (left); zoomed in loop route driven during the experiment (right). For details of the marked points with ‘×’, cf. Table 1. The other labels in round brackets denotes different paths along the experiment and are shown in next figure. The arrow indicates direction of travel during the experiment. Underlying map presentation taken from google maps

Figure 4. Images of different streets through which the kinematic experiment is conducted. Street view of simrockstraße (a), große barlinge (b), marienstraße (c) and schiffgraben (d). All the streets shown lie in the southern part of hannover and are indicated in Figure 3 with same corresponding labels

Table 2. Summary of marked points on the map

Each receiver was recording GPS, GLONASS and Galileo pseudo-range, carrier phase and Doppler observations and their corresponding signal strength. The sampling rate was set to 1 Hz for all the receivers and an elevation cut-off angle was not applied while recording the data. The geodetic-grade JAVAD GNSS receivers used in the experiment have the capability to track multiple frequencies (L1, L2 and L5) simultaneously, while the HS u-blox receivers can only track the L1 frequency of each GNSS system.