Use of Ultra-Low-Bandwidth Tracking Loops for Improved Dual-Frequency RTK Positioning

2.1 PLL Transfer Function

To understand the tracking performance of the ULB-PLL, it is necessary to derive the transfer function. Based on Figure 1, we denote the equation for the phase discriminator as follows:

$e_{\varphi }\left(z\right)=\varphi \left(z\right)-\hat{\varphi }\left(z\right)$1

For the carrier rate, we have the following:

$f_D\left(z\right)=a\left(z\right)+c_1{e}_{\varphi }\left(z\right)+g\left(z\right)/z$2

The feedback loop L1 is described by the following equation:

$g\left(z\right)=g\left(z\right)/z+c_2{e}_{\varphi }\left(z\right)$3

and the phase accumulation loop L 2 within the NCO is written as follows:

$\hat{\varphi }\left(z\right)=\hat{\varphi }\left(z\right)/z+T_{coh}{f}_D\left(z\right)$4

These equations can be rearranged to obtain the following closed form:

$\hat{\varphi }\left(z\right)=z{T}_{coh}\frac{a\left(z\right)(z-1)+(c_2-c_1+c_{1}z)\varphi \left(z\right)}{1+z\left(T_{coh}\right(c_2-c_1)-2)+z^2(1+c_1{T}_{coh})}$5

From this, a transfer function of the aided second-order PLL can be introduced as follows:

$\hat{\varphi }\left(z\right)=\left(H_{PLL}\right(z),H_{aiding}(z\left)\right)\cdot \left(\begin{array}{c}\varphi \left(z\right)\\ a\left(z\right)\end{array}\right)$6

with:

$H_{PLL}\left(z\right)=\frac{z{T}_{coh}(c_2-c_1+c_{1}z)}{1+z\left(T_{coh}\right(c_2-c_1)-2)+Z^2(1+c_1{T}_{coh})}$7

and:

$H_{aiding}\left(z\right)=\frac{z{T}_{coh}(z-1)}{1+z\left(T_{coh}\right(c_2-c_1)-2)+z^2(1+c_1{T}_{coh})}$8

To compute the equivalent noise bandwidth B_PLL, we can assume a signal without user dynamics, and a(z) can be assumed to be zero. By doing this, we verified that the coefficients for realizing an optimally damped loop are $c_1=1.89\sqrt{2}{B}_{PLL}$ and $c_2={\left(1.89B_{PLL}\right)}^2{T}_{coh}$, following the well-known GNSS tracking loop discussion summarized in Table 15.4 from Morton (2021).

To achieve a stable tracking loop, a stability condition must be fulfilled; this condition is commonly expressed by the requirement that the product B_PLL with update interval T_coh be much smaller than unity (Kaplan & Hegarty, 2017; Song et al., 2023). This analysis becomes complex, considering that a moving average filter is introduced for H_SAP in Section 2.5 and that the system has two inputs (true carrier phase ϕ(z), aiding signal a(z)). A stability analysis is left for future work, but it is noted that in this work and in the work by Dampf et al. (2025), the loop update rate T_coh and the beam-forming interval are short enough that their product with B_PLL is much less than unity.

2.2 Doppler Aiding

The aiding signal exploits information from an inertial sensor to estimate the LOS dynamics and is injected as a carrier rate signal in Figure 1. Thus, this method is commonly called Doppler aiding (DA). This signal is quantified as the aiding Doppler frequency a and consists of two components:

$a=f_{D,INS}+f_{D,CLL}$9

The first signal f_D,INS captures the LOS dynamics governed by the satellite v^S and user velocity v_INS as follows:

$f_{D,INS}=\frac{1}{\lambda }(v_{INS}-v^s)\cdot 1^e$10

where λ is the carrier wavelength and 1^e is the unit vector between the receiver and satellite expressed in the Earth-centered Earth-fixed frame (indicated by the superscript e). The second component, i.e., f_D,CLL, is related to the estimated receiver clock drift and is described in Section 2.3.

2.3 Clock-Locked Loop

As indicated in Equation (9) and pointed out a long time ago by Zhodzishsky et al. (1998), the receiver clock and multipath components remain in the signal after the LOS dynamics have been removed and need to be tracked by the PLL/DLL. The CLL realizes a centralized tracking loop for all tracking channels to track the receiver clock component. The CLL is realized as a second-order tracking loop, similar to a PLL/DLL (as shown in Figure 1). There are several fundamental approaches for estimating the receiver clock from the discriminator values, including (a) averaging all channel discriminator values at a specific time instant (Bochkati et al., 2023; Feng et al., 2022), (b) employing a CLL to obtain filtered estimates from the noisy discriminator values (Bochkati et al., 2023), and (c) using a Kalman filter (KF) for optimal estimates when a clock model is available (Lashley & Bevly, 2013b). This work adopts approach (b), utilizing a CLL with additional outlier checks, as shown in Figure 2. The implemented algorithm follows a two-step approach: in the first step, the relative clock error and clock drift are estimated using the CLL; in the second step, the clock drift is incorporated into the aided DLL/PLL. Each processing step, denoted as a “run,” postprocesses either a sequence of data or the entire data set. The first and second runs are indicated by black and red lines, respectively, in Figure 1 and Figure 2. In the first run, the input from the CLL, i.e., the receiver clock drift shown as f_D,CLL in red in Figure 1, is zero but becomes estimated as ${\hat{\dot{\varphi }}}_{CLL,i_s}$, as shown in Figure 2. In the second run, the preprocessed clock drift $f_{D,CLL}={\hat{\dot{\varphi }}}_{CLL,i_s}$ is applied in Figure 1.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Block diagram of the clock error tracking concept with time synchronization and outlier checks, before the data are filtered with a second-order CLL

The CLL algorithm takes as input the carrier-phase discriminator values from all tracking channels. These values are expressed in meters, making them independent of the carrier frequency. The carrier-phase error, derived from the carrier discriminator, is then provided as input to the CLL as follows:

$e_{\varphi }\left(i_j\right)=\varphi \left(i_j\right)-\hat{\varphi }\left(i_j\right)$11

where i_j denotes the time index for tracking channel j.e_ϕ,i,j = e_ϕ(i_j) is the short version, as shown in Figure 2, and the schema applies for all other variables. ϕ(i_j) is the raw carrier phase, and $\hat{\varphi }\left(i_j\right)$ is the estimated carrier phase. All three terms are shown in the z-domain by e_ϕ(z),ϕ(z), and $\hat{\varphi }\left(z\right)$ in Figure 1. The carrier discriminator values rely on the individual integrate-and-dump cycles of the respective channel correlator; thus, they must first be synchronized to a common time. The time sync block in Figure 2 implements the time synchronization step. The common interval over all channels is aligned to the first occurring epoch t_rx,j=1,i=1 of the first active tracking channel in the receiver time scale aligned to full milliseconds and is increased by the lowest coherent integration time T_ch,min in full milliseconds of all channels as follows:

$t_s=t_{rx,j=1,i=1}+i_s{T}_{coh,\min }$12

Hereby, the individual discriminator values at time index i_j are linearly interpolated to the common synchronized time index i_s by using the rate estimates of the respective tracking loop, as follows:

$e_{\varphi ,j}\left(i_s\right)=L\left(e_{\varphi }\right(i_j),i_s)$13

where L(…) denotes a linear interpolation function to the common time index i_s. Here, e_ϕ,is,j = e_ϕ,j(i_s), as shown in Figure 2. The next step is a carrier-to-noise power density (C/N₀) and elevation cut-off filter, shown as separate blocks in Figure 2, which filters the observations as follows:

$e_{\varphi ,k}\left(i_s\right)\subseteq e_{\varphi ,j}\left(i_s\right)\forall jwhere(s_{C/N_{0,}j}>{\tau }_{C/N_0})and({\theta }_j>{\tau }_{\theta })$14

Here, e_ϕ,k(i_s), k ∈ [1, K] is a filtered subset of e_ϕ,j(i_s), j ∈ [1, J] with $\kappa =\{j_{k=1},\dots ,j_{k=K},\}$ as a one-dimensional index look-up table of size (Kx1) to map the active channel index k to the global channel index such that $j=\kappa \left(k\right)$. Furthermore, s_C/N0,j is the signal-to-noise ratio, and θ_j is the satellite elevation angle of channel j, with τ_C/N0 and τ_θ as the respective filter thresholds in dBHz and radians, respectively. These two parameters reflect that the performance of the proposed CLL method is primarily determined by the dilution of precision and the availability of favorable satellite geometry. In addition, we employed the squaring C/N₀ estimator as described by Pany and Eissfeller (2006). In this method, C/N₀ is estimated blockwise for each preprocessed sufficient statistical data dump, on which the receiver operates. The selected block length for the processing setup was 500 ms. Because this estimated parameter is not further utilized, it does not have an impact on the results presented here. The outlier detection and phase realignment block filters out unreliable carrier-phase observations, ensuring that the CLL estimation algorithm is based on robust data. The following equation calculates the carrier-phase difference between the discriminator value e_ϕ,k(i_s) (expressed in meters) and the CLL phase:

$\delta \varphi (i_s,k)=2\pi \left(e_{\varphi ,k}\right(i_s)-{\hat{\varphi }}_{CLL}(i_s\left)\right)/{\lambda }_k$15

where δϕ(i_s, k) is expressed in radians, ${\hat{\varphi }}_{CLL}\left(i_s\right)$ is the estimated relative clock error from the CLL in meters, and λ_k is the wavelength in meters of the signal tracked by channel k.

For each channel, the CLL maintains a reference phase offset, which keeps the carrier-phase difference of the tracking channel and the estimated CLL near zero to avoid cycle-slip issues within the CLL. The initial value for the reference phase offset of each tracking channel is set as follows:

$\delta {\varphi }_{offset,k}=\delta \varphi (i_s,k)$16

and should be set after the PLL has converged. This step can be achieved by a delayed initialization (delay time of approximately 1/B_PLL seconds, where B_PLL is the PLL loop bandwidth in Hz).

For each channel k, the CLL continuously calculates the aligned phase difference:

$\delta {\varphi }_{aligned,i_s,k}=\arg \{e^{i\delta \varphi (i_s,k)}\cdot e^{-\delta {\varphi }_{offset,k}}\}$17

where arg{} delivers the angle of the complex number in radians and i is the imaginary part. The CLL monitors δϕ_aligned,k for phase outliers and resets δϕ_offset,k for the current i_s by applying Equation (16), if the following condition occurs:

$\left|\delta {\varphi }_{aligned,i_s,k}\right|\ge {\tau }_{{\varphi }_{outlier}}$18

where τ_ϕoutller = π/4 is an empirically evaluated and defined threshold for resetting the phase alignment value.

The relative clock error is estimated from the aligned discriminator values by the mean over all signals as follows:

${\overline{\delta \varphi }}_{CLL,i_s}=\frac{1}{2\pi K}\sum _{k=1}^K\delta {\varphi }_{aligned,i_s,k}\cdot {\lambda }_k$19

It should be noted that it is essential to express Equation (19) in meters, as we are combining phase observations from different frequencies.

In the next step, the noisy relative clock error estimate ${\overline{\delta \varphi }}_{CLL,i_s}$ is filtered with a well-known second-order loop, similar to that described in Equations (1)–(5) and further outlined by Dampf et al. (2022) and Teunissen and Montenbruck (2017), with a loop bandwidth of B_CLL = 40 Hz. The obtained output of the loop filter is the clock drift ${\hat{\varphi }}_{CLL,i_s}$, which is integrated using a NCO to the relative clock error ${\hat{\varphi }}_{CLL,i_s}$. As shown in Figure 2, the clock drift is, finally, forwarded to the individual tracking channels for aiding purposes, and the relative clock error is used for the outlier detection and phase realignment step in the CLL itself.

The clock estimation process is based on carrier-phase observations; thus, the estimation algorithm requires proper filtering such that only high-quality observations are used. Furthermore, it can be expected that in very challenging environments without any meaningful carrier-phase observations, the algorithm will not work well. In contrast, the high number of currently available signals from different GNSS systems and frequencies is very beneficial, because there is a higher probability that at least some good observations remain for the estimation process. It is worth noting that there is still potential to optimize the realized CLL algorithm, e.g., by using high-quality clocks with known behavior, but this will not be further discussed in this work.

2.4 Data Wipe-Off and DLL Aiding

The CLL algorithm discussed in Section 2.3 requires the use of a pilot signal. While Galileo provides pilot signals on all frequencies and GPS provides a pilot signal on the L5 frequency, the data bits must be removed from the GPS L1 C/A code signals. Thus, a data wipe-off was performed on the L1 C/A signal by using GPS C/A data bits from GFZ (Beyerle et al., 2008) based on the description by Kaplan and Hegarty (2017). The GFZ navigation data messages are open to the scientific community and cover 24 h of navigation messages for all available GPS satellites since June 18, 2007.

To improve DLL tracking loops, external aiding is nearly as vital as PLL aiding. Fortunately, the same aiding signal (Equation (9)) can be employed by scaling it with the ratio of the nominal code rate to the nominal radio frequency (RF) center frequency.

2.5 Thermal Noise and Carrier Multipath Mitigation

Ward (2017) recalled in Equation (8.70) that the thermal noise tracking error is as follows:

${\sigma }_{PLL}^2=\frac{B_{PLL}}{C/N_0}(1+\frac{1}{2C/N_0{T}_{coh}})$20

It is well known that the presence of the aiding signal a allows the PLL bandwidth to be reduced as LOS dynamics are removed from the burden of the tracking loop. Furthermore, it has been observed that the aiding allows one to increase the coherent integration time, as the predicted carrier phase $\hat{\varphi }$ controls the carrier replica generation to obtain the prompt correlator values P.

Assuming precise aiding, multiple prompt correlator values will share the same phase error, and thus, K correlator values can be accumulated by, e.g., a moving average (MA) described in the z-domain as follows:

$H_{SAP;MA}\left(z\right)=\frac{1}{K}\sum _{k=1}^K{z}^{-k}$21

As we want to demonstrate in future publications, this simple moving averaging filter together with the aiding realizes an SAP (Dampf et al., 2025), with the resulting beam pattern directed toward the satellite. In this work, we only want to note that the effective coherent integration time changes as $T_{coh}\to K{T}_{coh}$, further minimizing the thermal noise error (Equation (20)).

As nicely shown in Figure 22.7 by McGraw et al. (2021), the carrier multipath is described by a phasor diagram. The LOS and multipath signals add up in the complex plane, with the relative phase determined by the additional delay of the multipath signal with respect to the LOS signal. Assuming a constant multipath delay change of v_MP in meters per second, we obtain an oscillation of the true carrier phase around the value governed by the LOS signal, with a frequency of v_MP/λ, which can be derived from the Doppler shift formula f_D = –v/λ. Because the aiding signal does not contain this oscillation, a filter effect is applied, thereby reducing the multipath contribution to the final carrier-phase estimate $\hat{\varphi }$. The resulting carrier-phase multipath error amplitude as a function of the frequency v_MP/λ is shown in Figure 3. For this figure, typical parameters values are used: T_coh = 4 ms, B_PLL = 0.2 Hz, and K = 125.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Relative amplitude of the carrier-phase multipath error, assuming a constant multipath delay increase of v_MP for an aided second-order PLL, with and without the application of a moving averaging for the prompt correlator values

The carrier-phase multipath error is expressed as a relative amplitude. A relative amplitude of 1 indicates that the full multipath error of Figure 22.7 in the chapter written by McGraw et al. (2021) is contained in the RINEX output defined by $\hat{\varphi }$. Figure 3 clearly shows that the ULB-PLL mitigates most of the dynamic multipath owing to its low bandwidth. This mitigation is further enhanced by applying the moving averaging filter, i.e., SAP. The multipath delay change v_Mp depends on the radial velocity of the satellite (e.g., geostationary satellites or satellites at the zenith will cause lower values) and on the user velocity (e.g., a static user will again see lower values for v_MP compared with a user moving at high speed); however, it should be noted that upcoming navigation satellite systems in the low Earth orbit show much higher LOS dynamics, thus more fully exploiting the mitigation potential of the ULB-PLL.

4.1 Description

This data set is based on error-free simulated IMU observations, i.e., accelerations and rotation rates, where the ground-truth trajectory (see Figure 5) is also generated by the same tool. Here, the ground-truth, i.e., PVT, is given at sampling rates of 200 Hz, which is sufficient for use as an external DA source at similar sampling as the coherent integration time used in the MSRx (4 ms). Therefore, there is no need for the time-consuming fusion of GNSS/INS observations to obtain a high-quality trajectory for the DA of the tracking loops, as required by our proposed technique. By utilizing the available transceiver functionality of our software receiver (Pany et al., 2019), intermediate-frequency (IF) samples for GPS L1 C/A, L5Q and Galileo E1C, E5aQ have been simulated along the given reference trajectory. To achieve a better understanding and tuning of the CLL tracking techniques, as previously introduced by Bochkati et al. (2023), we introduced in the simulation both clock bias and drift as sine and cosine oscillations with amplitudes of 10^–9 s and 10^–9 s/s, respectively, at a frequency of 0.2 Hz. To enable RTK positioning, a simulation of a GNSS static reference station was conducted by applying the same settings, i.e., same Extended Standard Product 3 orbits and atmospheric corrections, as the rover. Here, the main advantage is the error-free time synchronization of the DA information and the absence of lever-arm uncertainty between the GNSS antenna phase center and IMU origins. For this reason, an overall 1- to 2-cm RTK position accuracy with 100% fixing of the carrier-phase ambiguities is expected, and hence, any degradation of the performance can be potentially linked to the quality of the implemented aiding techniques. Consequently, the hardware setup, as a potential error source, can be excluded from the error budget of the MSRx processing. For this contribution, we will consider 100 s of the simulation, which is sufficient to validate our concept.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Overview of the simulated trajectory, with the starting location at the University of the Bundeswehr Munich. (a) 2D view of the simulated trajectory, (b) Simulated velocity components (north v_N, east v_E, and down v_D), expressed in the local navigation frame.

4.2 Evaluation of the ULB-PLL With Optimal Signal Conditions

For the sake of clarity, we selected one GPS satellite with medium elevation to analyze the ULB-PLL tracking parameters, which are used through the entire processing tool chain of the MSRx. The tracking performance is depicted in Figure 6, with further details in the figure caption. For verification reasons, the code-minus-carrier (CMC) parameter is also shown in order to identify erroneous behavior that could originate from either implementation issues or, later, the quality of the tracked satellite signal. In the case of this selected satellite, the CMC exhibits the expected behavior, where, for example, no cycle slip or drift is visible. This result indicates that all aiding information has been properly applied to the tracking loops. The noisy behavior observed from 0 s to approximately 5 s is due to the transition from standard PLL tracking to the aided counterpart, i.e., ULB-PLL. This transition behavior is also well known from VT implementations.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

DLL/PLL tracking loop performance for the simulated data set for the GPS satellite G06

The left side shows the internal code and carrier rate based on DLL/PLL NCO readings and the estimated CLL clock drift. On the right side, the total (=internal + aiding) carrier pseudorange rate, the external DA signal, and the CMC are depicted.

A mathematical description of the CLL has been given in Section 2.3, where the two-stage CLL estimation is applied to extract the oscillator drift, which is shared by all tracking channels. Owing to the fact that the CLL is able to capture the simulated oscillator dynamics with the correct amplitude and periodicity in the first run (see Figure 7(a) or the lower left part of Figure 6) and later extract it from the tracking signals (Figure 7), the functionality and performance of the CLL has been successfully validated. Similarly, the impact of the two-stage CLL estimation can be translated to the discriminator level, as shown in the heat-map representation in Figures 7(c) and (d) after the first and second run, respectively. The heat-map is a colored histogram of the phase discriminator values from all signals.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

(Left) Estimated clock bias, clock drift, and carrier discriminator heat-map (for L1/E1 and L5/E5a) after the first run of the CLL, which aims to estimate the clock error; (right) the second/main run, where the sine wave of the simulated clock error has been almost entirely eliminated, resulting in discriminator values that show only the expected noise level of the carrier phase. (a) Estimated clock and clock rate (gray lines) of the first run (with red and blue curves representing the L1-CLL and L5-CLL discriminator values, respectively), (b) Estimated clock and clock rate of the second run, (c) Phase discriminator heat-map after the first CLL run for all satellites with L1 (upper plot) and L5 signals (lower plot), where oscillations of the clock error are clearly evident, (d) Phase discriminator heat-map after the second CLL run, i.e., after the estimated clock error has been applied to the tracking loops.

5.1 Description

The TEX-CUP data set (Narula et al., 2020) is a public benchmark data set collected in the dense urban center of the city of Austin, TX, which is dedicated for the evaluation of multi-sensor GNSS-based urban positioning with various challenging surrounding conditions. In contrast to other open-source data sets, this data set provides raw signal samples recorded after the analog-digital conversion step. These wideband IF GNSS sample data are recorded together with tightly synchronized raw IMU data and stereoscopic camera data.

For this publication, we make use of the raw data provided by the following hardware components:

• NTLab B1065U1-12-X RF front-end, configured to capture L1/L2/L5 signals from one GNSS antenna, which is also provided with a Bliley LP-62 low-power 10-MHz oven-controlled crystal oscillator external reference.

• Triple-frequency (L1/L2/L5) high-performance GNSS patch antenna from Antcom (NGS code: ACCG8ANT 3A4TB1).

• High-performance industrial-grade MEMS IMU from LORD MicroStrain, 3DM-GX5-25 attitude heading reference system.

• Two Septentrio AsteRx4 receivers, which are multi-frequency, multiconstellation dual-antenna receivers. One receiver was used as a reference rover, while the second served as a reference station to allow differential correction.

5.2 Raw Data Preprocessing

To make the raw data of these sensors compatible with the MSRx, a preprocessing step is required. First, the provided IF samples from the NTLab front-end were decoded based on the “ION GNSS software-defined radio metadata standard,” as described by Gunawardena et al. (2021), and then converted to sufficient statistics (Bochkati et al., 2022) by means of the MuSNAT software (Pany et al., 2019). This procedure is also summarized in the left block of Figure 4. A new GNSS/INS reference trajectory was generated based on the tight coupling of the Septentrio receiver observations and the LORD MicroStrain IMU signals, where the satellite observations from the reference station have been used to perform RTK positioning. The resulting PVT information was made available at a sample rate of 100 Hz to provide a high-quality DA signal for the MSRx.

Owing to the high computational load, only the first 400 s of the TEX-CUP trajectory were selected to evaluate our proposed technique. This trajectory snapshot is depicted in Figure 8 and highlighted by an orange rectangle. The beginning and end of the vehicle track are indicated by a green and cyan circle, respectively. Here is a brief description of the evaluated segment of the trajectory: this segment begins in an open-sky environment, where the measurement campaign commences. The vehicle then moves toward the city, as indicated by the blue arrows, encountering urban canyon conditions for approximately 2 min (see Figure 9). During this time, the vehicle passes under a small bridge and stops at a traffic light twice (as indicated by the red rectangle near the end of the trajectory). Significant signal degradation was anticipated at this location, which could pose a challenge for our ULB-PLL tracking technique. Furthermore, the satellite availability diagram shows good coverage of both GPS and Galileo satellites during data collection. It is important to note that the GPS L5 observations were not tracked by the Septentrio receiver at the reference station. An attempt to recover the L5 signal stream from the additional available native IF samples of the reference station was not successful. For this reason, the RTK evaluation of this data set is based on only the GPS L1 C/A and Galileo E1C/E5aQ signals. GPS L5 signals are, however, tracked with the ULB-PLL.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

(Left) Overview of the selected trajectory portion from the TEX-CUP data set; (center) skyplot of the available GPS and Galileo satellites; (right) Google street view imagery from one of the challenging scenarios selected in the TEX-CUP data set, marked in the map by a small red square

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

The multi-correlator plot for the GPS L5 frequency reveals two notable multipath signatures, highlighted within red ellipsoids. (a) Multi-correlator plot for GPS G04, (b) Zoomed-in image of Figure 9(a).

The first significant multipath occurrence is observed at 70–90 s, marked by magenta dots, whereas the second occurrence is noted at 210–230 s, indicated by cyan dots.

5.3 Performance of the ULB-PLL With Real Satellite Signals

To showcase the benefits of our advanced tracking technique, even under challenging signal conditions, the following evaluation strategy was applied. A reference/baseline MSRx configuration that applies standard PLL tracking was set up, where the output is compared with a second configuration with ULB-PLL, i.e., after DA and CLL have been applied. The key performance parameters for this comparison are the code- and carrier-phase residuals, which are estimated by RTKLIB. In other words, this software estimates RTK float residuals for both the baseline and ULB-PLL configurations, separately. Based on the simulation results, it is expected that the code/carrier residuals with the ULB-PLL will be considerably smaller than those of the baseline. The second key performance indicator is the deviation of the RTKLIB-estimated float position, for both configurations, compared with the TC-GNSS/INS reference trajectory. Similar to the analysis of the tracking loop outcome shown in the simulation part, i.e., in Figures 6–7, two different GPS and Galileo satellites with similar elevations and azimuths were selected for the case of the TEX-CUP data set. According to the skyplot in Figure 8, G07 (elevation ≈ 60°, azimuth ≈ 120°) and E11 (elevation ≈ 50°, azimuth ≈ 145°) coincide with this selection. The corresponding tracking loop parameters processed by the MSRx are depicted in Figure 10 and Figure 11, respectively. In both figures, the code rates (upper left plots) contain considerable spikes with different amplitudes, which can also be found in the CLL clock drift (lower left plots). This result is due to the fact that the estimated receiver clock error impacts all tracking channels. Additionally, for both selected satellites, the external DA is similar to the carrier rate, albeit smoother, which is also expected, because the total carrier rate also contains the jitter from the carrier NCO and the discriminator. As pointed out before, the correctness of applying both CLL and DA to realize the ULB-PLL technique can be roughly summarized in the CMC parameter. Ignoring the offset from zero, the CMC curve exhibits the expected behavior. The small peaks (with an amplitude of approximately 2 m) that are visible correspond to those present in both the code rates and CLL signal. The latter parameter is analyzed in Figure 13, where the outcomes of both estimation stages (see Section 2.3) are depicted. As shown in this figure, only a few satellites were engaged in the CLL, owing to the 50° elevation mask filter employed to reject noisy signals of the lower satellites from the CLL estimation process. Both jumps at approximately 70 and 210 s are shared by all used satellites, as shown in Figure 13(b). However, the CLL was able to detect these jumps correctly, as shown by the gray curve. Figure 12 combines the MSRx tracking parameter from ULB-PLL and the baseline configuration for GPS G04, where the improvement achieved by our tracking technique is evident. The combined impact from both CLL and DA is more visible in the code rate (upper left plot) and carrier rate (upper right plot), especially after 329 s, where the orange curve from the ULB-PLL smoothly crosses the noisy signal from the baseline. This behavior is also reflected in the CMC plot, where the amplitudes of many spikes have been considerably reduced.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 10

DLL/PLL tracking loop parameters from the TEX-CUP data set for a selected GPS satellite G07

See Figure 6 for an explanation.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 11

DLL/PLL tracking loop parameters from the TEX-CUP data set for a selected Galileo satellite E11

See Figure 6 for an explanation.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 12

DLL/PLL tracking loop parameters for G04, which show the improvements achieved by ULB-PLL (orange curves) over the baseline configuration

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 13

Two-stage CLL outcome showing both the clock bias $\hat{d}t$ and drift $\hat{\dot{d}}t$ estimated for the trajectory snapshot taken from the TEX-CUP data set; Both runs use a second-order CLL with a 40 Hz loop bandwidth (Bw). (a) First-run CLL, (b) Second-run CLL.

As highlighted in Figure 4, after the code- and carrier-phase information have been generated for both MSRx tracking configurations, two RINEX files can be obtained. Later, an RTK solution is computed with the RTKLIB software (version demo5 b34h) (rtklibexplorer, 2020), where only a float ambiguity solution is intended. The reason for selecting a float solution is that we want to obtain statistical measures of the improvements, which is difficult to achieve with fixed solutions having a binary nature. These measures can be obtained by setting a higher threshold for the ambiguity resolution, e.g., 999. From the RTKLIB output files, we used the *.stat-files, which contain both code- and carrier-phase residuals. In Figures 14(a) and 14(b), the estimated residuals from the two configurations are depicted. The different colors refer to the tracked satellites; however, for the sake of a better visualization, the pseudorandom noise identifier has been disabled from the legend. A noticeable improvement in the ULB-PLL residuals is evident at first glance. In addition to a reduced number of outliers, the residuals also appear smoother. A zoom into these signals, as shown in Figures 14(c) and 14(d), confirms the previous statement, where the bias of the residuals almost disappears and the random behavior for the baseline configuration follows an almost linear flat trend. If we translate this reduction in the residual level in numbers, we can see that, for example, the RMS error has been reduced from 2.23 m to 2.06 m for the code observations and from 0.12 cycles to 0.10 cycles for the carrier-phase observations. More statistics for these parameters are listed in Table 3. Because it is not straightforward to link the reported improvement in the residual domain to the position domain, the cumulative distribution function (CDF) and east-north-up (ENU) absolute position errors compared with the ground-truth trajectory have been computed. For this task, the same RTKLIB residual output file (*.stat) was used to extract the estimated position and compare it with the high-quality reference trajectory. The obtained results of the ENU position error and the corresponding CDF are shown in Figure 15(a) and Figure 15(b), respectively. According to the computed RMS value (shown in the legend), the ULB-PLL method achieved 0.28 cm in the east and 0.51 cm in the north component, corresponding to 39 cm and 22 cm, respectively, better than the benchmark configuration. Despite the offset visible in the up direction for both configurations, the ULB-PLL shows a 0.5-m improvement in the RMS value. This difference can be classified according to the probability by using the CDF representation (see Figure 15(b)). For this reason, we computed the error probability at 90%, 95%, and 99%, which covers the two-dimensional (2D) horizontal domain and 3D domain while accounting for the height error probability (Table 2). This table clearly demonstrates that the location reported by the ULB-PLL is less than 0.975 m away from the true location for 95% of the time, whereas the baseline configuration shows a distance of 1.345 m for the same probability distribution. This result reflects an improvement of approximately 70%. Similar behavior can also be observed in the 3D domain. However, the most impressive improvement is related to the 99% threshold, where our proposed technique achieved 1.140 m vs. 2.477 m for the 2D CDF error probability, which demonstrates an integrity improvement due to the use of the ULB-PLL.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 14

RTK float residuals estimated by RTKLIB with the MSRx baseline and ULB-PLL configurations. (a) Carrier- and code-phase residuals resulting from the baseline configuration, (b) Carrier- and code-phase residuals resulting from the ULB-PLL configuration, (c) Zoomedin plot of Figure 14(a) between 130 and 150 s, (d) Zoomed-in plot of Figure 14(b) between 130 and 150 s.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 15

(Left) Position error computed from the float solution obtained by RTKLIB (L1/E1/E5a) and the employed reference trajectory (TC Septentrio receiver/MicroStrain-IMU), with the positioning error expressed in the local ENU frame; (right) 2D/3D CDF of both the baseline and ULB-PLL MSRx configurations computed relative to the reference trajectory of the TEX-CUP data set. (a) ENU position error, (b) 2D/3D CDF.

5.4 Quality RTK Assessment of Carrier-Phase Observations

In addition to the evaluation of the RTKLIB float residuals, we processed the resulting RINEX data with commercially available software from NovAtel (NovAtel Inc., 2018) using their advanced RTK engine. Again, we used the available forward processing mode without the integration of any IMU data into this processing step. This mode refers to a KF-based forward estimation, i.e., without backward smoothing of the observation. The nearby reference station provided only GPS L1 and Galileo E1/E5a observations, with no L5 data. In addition to the baseline and ULB-PLL configurations, this comparison also includes a RINEX observation file from the reference GNSS receiver (Septentrio), which is a part of the reference trajectory (see Section 5.2). Therefore, the results from the Septentrio receiver can be viewed as a performance threshold in this RTK evaluation.

After processing three appropriate RINEX files by applying identical settings, we chose two parameters to asses the impact of the ULB-PLL technique. The first parameter is the ambiguity status, as depicted in Figure 16. The second parameter is described by the ENU position error with respect to the available ground-truth trajectory, as visualized in Figure 17. Specifically, in the case of ULB-PLL tracking, all carrier-phase observations have been accepted in the ambiguity estimation procedure, i.e., without the rejection of any observations. This result can be seen in Figure 16(c), indicated by the continuous green-colored surface. In this case, the fixing of the carrier ambiguities was successful for the entire data set, corresponding to 100% ambiguity-fixing success rates. This result indicates a good quality of the ULB-PLL-generated carrier-phase information. In contrast, the observation of the geodetic-grade Septentrio receiver showed surprisingly degraded quality, because some epochs were ignored (areas with background color) by the Inertial Explorer (IE) software and a few epochs were declared as the float solution (indicated by red bars). Table 4 summarizes the content of Figure 16 in numbers, where the amelioration of GNSS observation by applying the ULB-PLL technique can be clearly seen. The satellite observations from the baseline configuration yield 92% ambiguity-fixing rates, while two epochs are indicated as float solutions and six epochs are ignored. Similarly, the Septentrio receiver achieves the same ambiguity-fixing rates as the baseline, while two more epochs are engaged in the float solution and only four epochs are rejected.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 16

Carrier-phase ambiguity resolution status computed by the RTK engine of the NovAtel IE. (a) Ambiguity status when using the RINEX observation file generated by the MSRx baseline configuration, (b) Ambiguity status when using the RINEX observation file available from the reference Septentrio GNSS receiver, (c) Ambiguity status when using the RINEX observation file generated by the MSRx ULB-PLL configuration.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 17

Positioning error of the baseline, ULB-PLL, and reference Septentrio GNSS receiver compared with the TC Septentrio/MicroStrain-IMU reference trajectory. (a) East position error, (b) North position error, (c) Up position error.

This remarkable improvement in the quality of the carrier-phase observation — and of the code observations — is reflected in the absolute position accuracy. Here, RTK-based position results from the IE RTK engine were compared with the ground-truth trajectory. The obtained position error is depicted in Figure 17, where, at first glance, the errors from our aiding techniques (brown curve) appear to be remarkably smaller. While both the baseline MSRx configuration and the Septentrio receiver show many outliers with an amplitude of up to 5 m, the ULB-PLL curve exhibits a smooth, robust course, with a maximum error of 18 cm and an RMS of 9.5 cm in the horizontal level. It must be mentioned that the y-axes in these plots have been limited to ±0.5 m in order to properly visualize the behavior of the curves. Therefore, a visualization of the outliers was excluded from these figures. Compared with the baseline configuration, the improvement through the ULB-PLL processing is up to three fold in the horizontal level, as shown in Table 5. In the same way, the RMS values derived from our methods are approximately 2 cm smaller. This parameter for the up-component has also been considerably improved, i.e., from approximately 0.5 m down to 5.7 cm.

A detailed analysis with this robust commercial GNSS positioning engine again shows that a substantial improvement is achieved by using the ULB-PLL compared with conventional DLL/PLL tracking. Furthermore, the ULB-PLL also outperforms the commercial rover receiver that is included in the TEX-CUP data sets.