Abstract
Precise measurements of signal phase are essential for Global Navigation Satellite System (GNSS) position estimates. However, propagation through the earth’s ionosphere imposes frequency-dependent phase errors. The frequency dependence is exploited to correct phase errors proportional to total electron content (TEC) divided by frequency. Scintillation causes additional stochastic errors, which can become the largest phase error contribution. Both geometry-free (GFCs) and ionosphere-free (IFCs) frequency combinations are subject to such uncorrectable but generally small phase scintillation errors. A recent published study compared GPS TEC estimates derived from nominally identical L1-L2 and L1-L5 GFCs. The TEC errors establish an upper bound on uncorrelated single-frequency scintillation error contributions. The measured TEC errors were verified with phase-screen simulations. However, a known but largely overlooked phase unwrapping error affects the extraction of signal phase from phase-screen complex signal realizations. This paper demonstrates a procedure for detecting and correcting phase-unwrapping errors and discusses their ramifications.
1 INTRODUCTION
Precise phase measurement is an integral part of GNSS signal processing. However, unlike signal intensity and delay, phase cannot be measured directly. Formally, phase is extracted from a complex product of two complex time series with nearly identical frequencies. In purely mathematical terms, phase is derived from a 2π-ambiguous nonlinear inverse tangent operation. For GNSS applications, smoothly varying phase reports are essential, whereby a considerable ongoing effort has been devoted to detecting and repairing cycle slips (Breitsch et al., 2020). There is a well-know algorithm that will recover unwrapped phase from a sequence of 2π ambiguous samples. The algorithm is not used directly for GNSS processing because it operates on an extended set of data samples. However, phase reconstruction is used to provide phase truth for evaluating cycle-slip recovery schemes and GNSS system performance under disturbed propagation conditions.
Theoretical scintillation models are being used to generate realizations of propagation disturbances that can be mapped onto signal realizations for direct performance analysis. In the process of investigating the phase structure generated by phase-screen simulations, we found a known but largely neglected problem with the phase unwrapping operation that reconstructs the wrapped inverse tangent phase. The unwrapping operation only produces jump-free phase reconstructions when the sample differences do not exceed π radians. We demonstrate a modification of the unwrapping procedure that will detect and remove residual phase jumps. The remainder of this paper will describe a modified phase reconstruction utility that will detect and correct the phase derived from complex signal realizations. Readers interested mainly in the phase unwrapping problem can skip to the Complex Signal Phase Extraction section.
Global Navigation System Satellite (GNSS) receivers measure signal intensity, phase, and delay as frequency-dependent time series. Irrespective of the complex modulation that defines the GNSS signals, measurements can be interpreted against the following parametric signal model, which was introduced in Rino et al. (2018):
1
2
3
The leading term, r(t) in Equations (2) and (3), represents the true propagation distance (range), fc represents the carrier frequency, and λc = c/fc is the carrier wavelength with c as the vacuum velocity of light. The ellipses indicate omitted bias and noise contributions. Delay and phase have been scaled to pseudo-range length units. If the total electron content is measured in units of 1016 electrons per square meter per meter of path length, the product cK = 40.3. The conversion from path-integrated electron content to phase in radians is
4
The following definition of the complex modulation identifies the stochastic TEC and scintillation components explicitly:
5
Equation (1) in Nguyen et al. (2019) also identifies the ionospheric error terms explicitly. In their notation, Lf = φ(t; fc)λc/(2π). Our TEC and scintillation terms when scaled to range units coincide respectively with and in their notation. In textbook versions of the signal model, phase scintillation is often lumped together with other error terms.1 The error terms include noise, interference, multipath, and system biases as enumerated in Nguyen et al. (2019). The measured signal phase also includes a 2π phase ambiguity. For this study, only the ionospheric terms will be considered.
Geometry-free combinations (GFCs) of the frequency-dependent phase follow from Equation (2) as
6
In the absence of scintillation, GFCs estimate the combined deterministic and stochastic TEC contributions to the signal phase. Although single-frequency phase scintillation contributions cannot be isolated, the scintillation-induced TEC errors can be isolated by subtracting nominally identical TEC estimates derived from GPS L1-L2 and L1-L5 GFCs. A study reported by Rino et al. (2019) found TEC scintillation errors to be less than one TEC unit when the S4 scintillation index at the L2 or L5 frequency was less than 0.5. This implies that for the weak-to-moderate scintillation conditions that apply to most GNSS observations, scintillation-induced TEC errors for ionospheric diagnostic measurements are negligibly small. The TEC error bounds were corroborated in Rino et al. (2019) with phase-screen simulations.
2 PHASE-SCREEN PROPAGATION MODELS
The time series that characterize the complex modulation imparted to GNSS signals are generated by translation of the propagation paths through structured regions. A complete calculation requires knowledge of the varying positions of the source and receiver, the location and extent of the structure, the structure drift, and the direction of the magnetic field. Fortunately, a comparatively small number of parameters capture the critical dependencies on structure, propagation distance, and time-to-space translation. Propagation effects are determined by the Fresnel scale
7
where rF is the effective propagation distance from the structured region to the receiver, and k = 2πfc/c is the propagation vector magnitude. Space-to-time conversion is achieved by replacing the Fresnel-scale with ρF/vveff where vveff is an effective scan velocity, which incorporates the dependence on the scan direction relative to the magnetic field direction. Representative values can be used without specifying the propagation geometry details, as described in Carrano and Rino (2016), Jiao et al. (2018), Xu et al. (2019), and Carrano et al. (2019). A phase-screen model that captures the signal structure is described in Rino et al. (2018).
The model is revisited here because, as will be shown in the next section, phase-screen sampling impacts phase reconstruction. The phenomenon is demonstrated with complex field realizations initiated with stochastic phase realizations. The phase structure is characterized by a spectral density function (SDF), which is the average intensity of a Fourier decomposition of the process. For a time series, the term power spectral density (PSD) usually replaces SDF, which is the more general term. A standard procedure for generating phase realizations with SDF Φm follows:
8
The random variable, ηm, is a zero-mean, unit-variance, uncorrelated (white noise) process. Formally, 
. Angle brackets denote ensemble averaging. With the SDF, Φm, specified, a complex field realization can be generated with two discrete Fourier transformation (DFT) operations as follows:
9
10
11
where 𝜇m = ρFqm is a scale free measure of the spatial wavenumber.
Carrano and Rino (2016) developed a two-dimensional scintillation theory for unbounded two-component powerlaw structures. Efficient algorithms are available for computing the intensity SDF, ΦI, the scintillation index, S4, and the complex field SDF, Φψ. Calculations are initiated by specifying the following scale-free form of the two-component phase SDF:
12
where 𝜇0 = ρFq0 and 
. The defining parameters are Cpp, p1, 𝜇0, and p2. However, it was found to be more convenient to use the universal strength parameter
13
which determines the severity of the scintillation as measured by the S4 index
14
For 𝑈 < 1, S4 ≪ 1. As 𝑈 increases, S4 approaches unity. As 𝑈 increases further, S4 may exceed unity or achieve a local maximum after which S4 recedes to a value approaching unity from above.
The theory also provides analytic representations of the phase structure function, D, the mutual coherence function, MCF, and the complex field SDF, Φψ:
15
16
17
The notation δy is used to avoid confusion with the sample interval Δy. The MCF and the complex field SDF are defined by the phase screen. However, they remain invariant as the complex field propagates away from the phase screen. This well-known property of freely propagating stochastic plane-wave fields follows from the diffraction operator in Equation (11), which modifies only the phase of the propagating plane waves that comprise the field.
As discussed in Carrano and Rino (2016), the software implementation of the theory generates normalized SDFs as functions of 𝜇. The 𝜇 range is sampled non-uniformly from 𝜇min to 𝜇max. For applications, a uniformly sampled 𝜇 range is generated by interpolation. Physical units are defined by the spatial extent 𝑌 and the sampling interval Δy or the temporal extent T and the sampling interval Δt. The corresponding frequency ranges are
18
19
The Fresnel scale factors ρF and ρF/vveff convert the physical frequency ranges qm and fm to the scale-free frequency range 𝜇m. Hereafter, only spatial units will be used with the understanding that conversion to time and temporal frequency is achieved by replacing ρF with ρF/vveff.
The components of the Fourier transform pair MCF(η) ⇔ Φψ(q) provide guidance for sampling. If they are well-contained within the respective intervals 𝑌 and Q, both functions can be computed accurately with an N-point DFT pair with sampling Δδη and Δq, corresponding to the Δy and Δq sampling of the complex field realizations. Furthermore, Fourier interpolation, which multiplies N by zero-extending the discrete Fourier transform accordingly, accurately represents the interpolated complex field samples. However, the intensity and phase of the complex signal are derived from nonlinear operations that impose more stringent sampling requirements. As described in the Appendix, the theoretical intensity SDF can be used to determine sampling requirements for intensity reconstruction. We have no analytic theory against which phase reconstruction can be verified, whereby we must rely on consistency checks as described in the next section.
2.1 Complex signal phase extraction
Extracting phase from a discrete complex series, ψn, starts with the wrapping operation
20
If there are phase changes such that |φn−1 −φn|≥ π, the transitions cannot be unambiguously resolved without further analysis. However, if |φn−1 −φn| < π, perfect reconstruction can be achieved with a well-known unwrapping algorithm. We will refer to any single-step phase change greater than one cycle (2π radians) as a phase jump to distinguish it from a cycle slip as defined in Nguyen et al. (2019). Following the development of Itoh’s algorithm as reviewed in chapter 1.3.1 of Ghiglia and Pritt (1998), we apply the wrapping operation to the differences of the wrapped phase. The result is
21
where k2 is the unwrapping error for the difference. If there are no phase jumps, Δk1 −k2 = 0. A phase-jump-free phase sequence is recovered by summing the unwrapped difference W{ΔW{ψn}} = Δφn. Although perfect reconstruction can be achieved if there are no phase jumps at the current sampling interval, the phase associated with critically sampled complex fields may not satisfy the jump-free criteria.
To introduce some notation, let N×nint represent the size of a zero-extended DFT. An inverse DFT applied to the zero-extended original DFT will repeat the original samples at intervals nint followed by nint −1 interpolated samples. Thus, nint = 1 implies no interpolation. Contributions at the Nyquist frequency must be zero, which imposes a strict periodicity constraint. However, phase-screen realizations automatically satisfy this constraint. To test the unwrapping operation at interpolation level nintN, we generate an interpolated sequence 
with (nint + 1)N samples and apply the following test to see if the common samples have changed:
22
The double equal sign operation generates ones where there is equality and zeros otherwise. If there are no disparities, the phase has been correctly unwrapped at level nint. The operation can be repeated with nint replaced by nint + 1 until Equation (22) is satisfied. With noise-free complex samples, the interpolation limit is determined by numerical resolution.
To illustrate unwrapping error detection and correction, a 200-km stochastic phase realization was generated with 8,192 samples, whereby the initiating phase is known. The turbulent strength parameter Cpp was increased until unwrapping errors were observed when the unit-intensity complex signal phase was unwrapped. The upper left frame of Figure 1 shows the phase realization (blue) with the unwrapped phase overlaid in (red). The disparity between the initiating phase and the unwrapped phase is a small fraction of the total phase. However, the absolute errors shown in the lower frame are significant. The nominal phase difference computed from Equation (21) is shown in the lower left frame (blue). The right frame shows the same results with nint = 2, which reproduces the original phase sequence without error. For a critically sampled complex field, which is intrinsically smooth at the finest sampling level, interpolation can be increased until perfect reconstruction of the original phase samples has been achieved.
Upper left frame shows initiating phase (blue) and unwrapped phase (red). The S4 L1 index at 350 km is indicated for reference. The lower left frame shows the unwrapped phase (blue) and the error (red). The right-hand frame shows the same displays with unwrapping phase correction
3 UNWRAPPING ERROR RAMIFICATIONS
We showed in the previous section that the phase of a complex series may be inadequately sampled for conventional phase unwrapping even if the complex signal itself is critically sampled. The complex signal amplitude has no direct impact on phase unwrapping. However, strong amplitude scintillation is associated with highly structured phase scintillation, which is more prone to unwrapping errors. Phase-screen realizations were generated as described in the previous section. The upper frame of Figure 2 shows the 𝑈 = 1 observation-plane intensity realization derived from the phase realization shown in Figure 1. The intensity is sampled at the minimum sampling that supports error-free phase unwrapping of the complex signal phase, namely 3N = 24386. The lower frame shows the uncorrected phase (red) and the corrected phase (cyan).
Upper frame shows the interpolated intensity realization (red). Lower frame shows initial (red) and corrected (cyan) phase. The starting distance is arbitrary
The upper frame of Figure 3 shows the SDF of the interpolated intensity realization. The pentagrams mark the Nyquist frequencies for nint =1, 2, and 3 sampling. The intensity content beyond the complex signal Nyquist frequency occupies only the range from the Nyquist frequency to twice the Nyquist frequency. However, that extension is not an extrapolation of the theoretic SDF. Increasing the sampling of the initiating phase screen is necessary to increase the intensity resolution. The lower frame in Figure 3 shows the SDF of the signal phase (blue) with the scintillation SDF overlaid (red).
Upper frame shows intensity SDF for 𝑈 == 1 realization. Lower frame shows the phase SDF (blue), the scintillation SDF (red), and the theoretical SDF for the initiating stochastic phase (green)
Figure 4 shows the occurrence of unwrapping errors plotted against 𝑈 at the L1 frequency. The lower L2 and L5 frequencies correspond to higher U values. The occurrence of unwrapping errors is confined to stronger scintillation conditions. Fewer than six iterations with typically no more than three iterations generated error-free phase unwrapping. Interpolation does not change the common intensity samples, but as noted above, increasing the spectral range coincident with the theory requires recalculation of the propagation effects.
Phase L1 and L5 unwrapping errors versus 𝑈 at L1
To demonstrate the effects of phase unwrapping errors, Figure 5 shows a recalculation of the TEC errors shown in Figure 5 from Rino et al. (2019) with unwrapping error detection and correction. The distance scale has been normalized to generate a zero-to-one scale for the S4 index at L1, which is overlaid in red. The black dots, which are computed as (nint − 1)/10 for display, show the interpolation required to correct unwrapping errors. The unwrapping error occurrence is consistent with the results shown in Figure 4. The first seven realizations required no interpolation. In effect, the occurrence of unwrapping errors coincides with significant TEC errors. The main conclusion, namely that significant TEC errors are negligibly small under weak to moderate scintillation conditions, is unchanged.
TEC errors versus 𝑈 derived from complex field realizations corresponding to L1 and L2 realizations at 350 km. The ordinate is scaled to the maximum distance to allow superposition of the S4 index in red. The black dots indicate the level of iteration required for error-free unwrapping
4 SUMMARY AND CONCLUSIONS
We showed that a standard procedure used to extract a phase from a complex series may not generate a phase series free from phase jumps, as illustrated in Figure 1. These residual phase jumps are caused by sampled-dependent phase changes greater than π radians. A critically sampled complex series supports interpolation that will ultimately eliminate phase changes greater than π radians. We demonstrated a successive interpolation procedure that detects and removes any residual phase jumps. Functionally, critical sampling means that the phase associated with the complex field has no discontinuities. This constraint does not preclude the very large phase changes associated with the so-called canonical scintillation deep fades. Noise-free phase-screen simulations generate very large phase changes, but when properly sampled, the phase variations are resolved (Carrano et al., 2013).
We showed also that the occurrence of phase jumps in phase-screen realizations is confined to large propagation disturbances. Revisiting simulations of TEC scintillation errors identified with three-frequency geometry-free GPS measurements (Rino et al., 2019), we found that phase jumps occurred under the same conditions that produced large scintillation-induced TEC errors. There was no impact on the critical finding, namely that scintillation-induced TEC errors are negligible for diagnostic measurements under weak to moderate scintillation conditions. In the same way, recent performance evaluations reported by Xu et al. (2019) are largely unaffected by uncorrected phase jumps because they are statistically similar to large but resolved phase changes. For future studies that use complex field simulations as truth, corrected phase extraction should be used.
Diagnostic phase measurements are most often generated by signal-processing procedures that track incremental phase or frequency changes. Cycle slips, which are highly processor and sample dependent, are attributed to the processor errors. The study by Nguyen et al. (2019) associated cycle slips with phase scintillation as distinct from the more smoothly varying correctable TEC contributions. Furthermore, signal processing operations are acutely sensitive to noise, which has not been considered. In a recent paper by Breitsch et al. (2020), a systematic procedure for identifying rapid local phase changes as cycle slips and demonstrating their dependence on noise and scintillation intensity was demonstrated. Simulated and measured cycle slips showed remarkably similar characteristics. Comparing corrected unwrapped phase sequences and measured phase sequences should provide insights into both the diffraction processes that generate rapid phase changes and signal processing operations that encounter the changes.
An updated MatLab scintillation model with examples of unwrapping error correction and sample adequacy checks can be downloaded from https://github.com/cusense-lab/gnss-scintillation-simulator
HOW TO CITE THIS ARTICLE
Rino C, Breitsch B, Morton Y, Xu D, Carrano C. GNSS signal phase, TEC, and phase unwrapping errors. NAVIGATION. 2020;67:865–873. https://doi.org/10.1002/navi.396
ACKNOWLEDGMENTS
The results in the paper can be reconstructed with the MatLab software provided. The reported research was supported in part by the Air Force Office of Scientific Research grant FA8650-14-D-1725.
APPENDIX
The algorithmic encapsulation of the two-dimensional phase-screen theory in Carrano and Rino (2016) generates normalized SDFs as functions of 𝜇 = ρFq. The 𝜇 range is sampled non-uniformly from 𝜇min to 𝜇max. For applications, the non-uniformly sampled range must be interpolated to
A1
This involves a choice of ρF, Δq, and N. If estimates Ŝ4 
are available, model parameters can be adjusted to match the theory. Figure A1 shows the non-uniformly sampled model output for the defining parameters in the title with a uniformly sampled overlay.
Model calculations of Intensity SDF versus 𝜇. Blue curve is non-uniformly sampled model output. Red curve is a uniformly sampled subset for a specified ρF
The red curve might represent the best fit range for the measured SDF and S4. However, for modeling purposes, we often have a specified N and Δq for which a range of scintillation strengths are of interest. One can then ask what range of ρF values are appropriate for realizations defined by fixed structure parameters. The effect of changing ρF with all the other parameters fixed is to slide the red curve in Figure A1 along the blue curve. It is intuitively clear that a good choice for ρF captures the largest values of the intensity SDF. This can be formalized by choosing ρF to minimize the difference between the sampled estimate of
A2
and the theoretical value. An algorithm for performing this computation is included in the revised phase-screen model. The red curve shown in Figure A1 was plotted using the algorithm calculation of ρF. The choice of the Fresnel scale is constrained by the phase-screen sampling. For example, the optimum sample range for smaller and smaller values of N will ultimately generate a realization too small to resolve the contributing structure.
Footnotes
Funding information
Air Force Office of Scientific Research, Grant/Award Number: FA8650-14-D-1725.
↵1 See Misra and Enge (2012), equations (5.46) and (5.47).
- Received February 12, 2020.
- Revision received July 11, 2020.
- Accepted July 16, 2020.
- Copyright © 2020 Institute of Navigation
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.