Research ArticleRegular Papers
Open Access

Efficient Signal Quality Monitoring of GNSS Signals Disturbed by Evil Waveforms

Fernando D. Nunes and Fernando M. G. Sousa
NAVIGATION: Journal of the Institute of Navigation December 2024, 71 (4) navi.668; DOI: https://doi.org/10.33012/navi.668

Abstract

Evil waveforms (EWFs) are anomalies in signals transmitted by a global navigation satellite system, provoked by electric malfunctions, that can significantly degrade the accuracy of the receiver’s position, velocity, and time solution. In this work, cross-correlation functions of a received signal disturbed by EWF distortion and the locally generated code signal are derived for threat models TM-A, TM-B, and TM-C, with expressions obtained for binary phase shift keying, binary offset carrier (BOC), and composite BOC pilot modulations. Closed-form expressions are offered in terms of sine integral functions for TM-A. A description is presented on how to use these results to assess the performance of EWF detectors through semi-analytic simulation, allowing the detectability and hazard regions to be determined in a computationally efficient way.

Keywords

1 INTRODUCTION

Evil waveforms (EWFs) are anomalies in global navigation satellite system (GNSS)-transmitted signals provoked by electric malfunctions that occur in the signal generators aboard GNSS space vehicles (Phelts, 2001; Selmi et al., 2020; Thevenon et al., 2014). The first reported EWF event occurred in 1993 with Global Positioning System (GPS) pseudorandom noise 19, inducing positioning errors of several meters (Shallberg et al., 2017). These threats manifest themselves in the form of an anomalous correlation peak (Phelts & Akos, 2006). Because the peak is distorted and displaced from the expected position, the processing of EWF-distorted signals by the receiver may lead to a significant loss of accuracy in the position, velocity, and timing solution, thus preventing its use in most applications (Gómez-Casco et al., 2022). Three types of EWFs are usually considered in the literature (ICAO, 2018; Pagot, 2016; Phelts, 2001): threat models TM-A, TM-B, and TM-C, which are associated with digital, analog, and digital plus analog distortion, respectively.

Signal quality monitoring (SQM) techniques have been developed to detect distortions in GNSS signals caused by multipath or spoofing attacks (Ali et al., 2014; Irsigler, 2008). By incorporating additional early-late correlators for monitoring purpose, distortions on the autocorrelation peak can be detected, and corresponding signals can be deweighted (or excluded) (Pirsiavash et al., 2017). Similarly, SQM algorithms may also be employed for detecting the presence of different types of EWF distortion (ICAO, 2018). The performance of an EWF detector may be measured in terms of the probability of not detecting a signal anomaly for a predefined probability of false alarm.

According to Pagot (2016) and Pagot et al. (2018), the main EWF detection algorithms are based on three tests: T1 (simple ratio), T2 (difference ratio), and T3 (sum ratio), which use a bank of early and late correlators. Basically, these tests measure the distortion of the code cross-correlation shape, such as the flatness of the main peak and the correlation asymmetry, for a certain value of C/N0 and produce a decision concerning the presence/absence of an EWF. The EWF detectability in the threat space and the corresponding hazard regions can then be determined by either using experimental data or resorting to simulations (Negrinho et al., 2021).

The occurrence of EWFs is rare (Shallberg et al., 2017); thus, it is difficult to obtain experimental data that allow one to reliably evaluate SQM algorithms in the presence of EWFs and additive white Gaussian noise (AWGN). Instead, Monte Carlo (MC) simulations may be advantageously applied to assess the performance of the algorithms. MC simulations require very little mathematical analysis and can be applied to almost any communication system, although at the expense of very long simulation run times, especially when the a priori probability of a given effect is very small. Consequently, semi-analytic (SA) simulations, in which analytical and simulation techniques are used together, are often proposed as a rapid alternative. These techniques can be applied because the Gaussian channel noise is not affected by the nonlinearities induced by the EWF. Although SA simulations require a higher level of analysis, the payoff is a significantly reduced run time (Tranter et al., 2004).

Following our previous work (Negrinho et al., 2021; Nunes & Sousa, 2022; Sousa et al., 2022), we derive mathematical expressions for signal cross-correlation functions that can be applied to different threat models and modulations in SA simulations. As processing AWGN with the correlators is a linear operation, the result is a set of statistically dependent, zero-mean Gaussian random variables that are added to the various correlator responses to the EWF-distorted waveforms, allowing us to obtain computationally efficient SQM algorithms.

The remainder of this paper is organized as follows. Sections 2 and 3 analytically characterize the multi-correlator outputs for TM-A, TM-B, and TM-C. Section 4 discusses the metrics and tests used by the GNSS receiver for EWF detection. Section 5 presents examples of SA simulations in terms of detectability and hazard regions, whereas Section 6 presents some final remarks.

2 THREAT MODEL A

TM-A was originally defined for coarse acquisition (C/A) GPS signals in which the positive chips have a falling edge that leads or lags relative to the correct end-time for that chip (ICAO, 2018). The threat model is associated with a failure in the navigation data unit inside the digital partition of the satellite.

2.1 Threat Model TM-A1

This threat model corresponds to the digital distortion 1 described by Pagot et al. (2016): a lead/lag on every falling transition after modulation by the code signal. For this model of distortion, only the lead/lag parameter d is required (d < 0 indicates a lead and d > 0 denotes a lag). Typical waveforms for TM-A1 with sine-phased binary offset carrier (BOCs)(m, m) signals are presented in Figure 1, where we assume that |d| < Tc (chip duration). This means that the 1 and −1 state durations are not equal, resulting in a distortion of the correlation function and inducing a bias in the pseudorange measurement equal to half the difference in the durations (Raimondi et al., 2012).

FIGURE 1
FIGURE 1

Typical waveforms for TM-A1 with BOCs(m, m) signals and |d| ≤ Tc

Let ρ(t) ≡ s(t – d) – s(t). The distorted signal is Embedded Image, where the error signal is as follows:

Embedded Image 1

with sign(x) equal to −1, 0, or +1 if x < 0, x = 0, or x > 0, respectively. We now define the cross-correlation operator between ϵ(t) and s(t) as follows:

Embedded Image 2

where T denotes the correlation duration. Thus, Rϵs (τ) = (1/2)[Rs (τ – d) – Rs (τ) + 〈|ρ(t)|s(t – τ)〉], where the autocorrelation of s(t) is Rs (τ) = 〈s(t)s(t – τ)〉 and 〈|ρ(t)|s(t – τ)〉 ≈ 0 (see Figure 1), leading to Rϵs (τ) = [Rs (τ – d) – Rs (τ)]/2. The cross-correlation between Embedded Image and s(t) is Embedded Image. The cross-spectrum corresponding to Embedded Image is then as follows:

Embedded Image 3

where Embedded Image is the power spectrum of s(t) and Embedded Image denotes the Fourier transform.

The receiver uses a bank of correlators, with the EWF-distorted signal being correlated with the set of local undistorted replicas. The processing is carried out in baseband using a low-pass filter with transfer function H(f), as shown in Figure 2. The generic correlator output is given by the following:

Embedded Image 4

where the filter output is as follows:

Embedded Image 5

FIGURE 2
FIGURE 2

Cross-correlation between the incoming signal y(t) and the local replica g(t)

Here, h(t) denotes the filter impulse response. The cross-correlation RA(·) accounts for the autocorrelation function of the (baseband) GNSS signal, the distortion due to the EWF, and the filtering effect in the receiver front-end. For simplicity, we consider the receiver’s locally generated signal to be the two-level subcarrier g(t) = s(t) for binary phase shift keying (BPSK) signals and g(t) = sBOC(t) for BOC and composite BOC pilot (CBOCpilot) signals. Thus, the correlation is carried out between the filtered distorted signal and the corresponding locally generated non-distorted and unfiltered version of the signal. However, in the case of transmission of the CBOCpilot signal, we perform the correlation with the two-level BOCs(1,1) carrier, with this simplified scheme resulting in a relatively small performance degradation. Thus, we have the following:

Embedded Image

or in the frequency domain:

Embedded Image 6

leading to:

Embedded Image 7

There are different methods for calculating RA(τ) using Equation (7):

  1. - Determine an analytic expression for the integral. The result will be exact, but the calculation may be complex, except for very simple transfer functions H(f) and ideal spreading codes with an infinite period.

  2. - Calculate the integral numerically. This method is approximate, with the accuracy depending on the length of the integration step.

  3. - Use the discrete Fourier transform (DFT) to compute Embedded Image from the cross-correlation samples Embedded Image, multiply Embedded Image by H(nΔf), and apply the inverse DFT (IDFT) to calculate RA(mΔτ). Here, Δτ and Δf are increments in the time and frequency domains, respectively. This method is also approximate, but its accuracy can be controlled by the size of the DFT sequence (Oppenheim & Schafer, 1975).

Consider now that g(t) = s(t) in Equation (7), which is valid for BPSK and BOCs modulations. Using Equation (3) yields the following:

Embedded Image 8

Henceforth, for simplicity, we assume that the low-pass filter in Figure 2 is ideal, with bandwidth B and zero group delay, such that H(f) = 1 for |f|< B and 0 otherwise. Let x = fTc (normalized frequency) and D = d/Tc (normalized EWF delay). We then obtain the following:

Embedded Image 9

Note that RA(–τ) = RA(τ) (even function) only if the EWF distortion is absent (D = 0).

2.1.1 BPSK(m) Signals

Let us define the triangular pulse ΛL(x – x0) as follows:

Embedded Image 10

with L denoting the pulse half-duration and x0 denoting the pulse center.

If the spreading code is an ideal long spreading code, the autocorrelation function of the BPSK(m) signals is Embedded Image (Betz, 2016), where the chip duration is Tc = Tc0/m, m = 1, 2,… and Tc0 = 10−3/1023 s is the chip duration of the GPS C/A code. The corresponding power spectrum is as follows:

Embedded Image 11

with sinc(x) = sin(πx)/(πx). Applying Equation (11) in Equation (9) leads to the following:

Embedded Image 12

The integrals in Equation (12) are of the following form:

Embedded Image 13

yielding:

Embedded Image 14

Appendix A explains how Is2c(a, b; L) can be calculated in terms of trigonometric and sine integral functions.

It can be shown that the peak of RA,BPSK(τ) occurs at τ = DTc/2 and that the function is symmetric relative to this point. Consequently, the receiver’s delay lock loop (DLL) will exhibit a range error approximately equal to cDTc/2, where c is the speed of light.

The cross-correlation RA,BPSK(τ) is displayed in Figure 3 for BTc = 1 and different values of D. The selected normalized bandwidth BTc = 1 corresponds, for instance, to BPSK(10) with a pre-correlation filter of bandwidth 2B = 20.46 MHz, which is the Galileo receiver reference bandwidth for E5a and E5b bands (European GNSS Agency, 2023).

FIGURE 3
FIGURE 3

Cross-correlation RA,BPSK(τ) for different lags D with BTc = 1; modulation: BPSK(m)

2.1.2 BOCs(m, m) Signals

The power spectrum of the BOCs(m, m) modulation is given by the following (Sousa & Nunes, 2013):

Embedded Image 15

Substituting Equation (15) in Equation (9) leads to the following:

Embedded Image 16

Figure 4 shows cross-correlation curves obtained from Equation (16) for different lags D and BTc = 10.

FIGURE 4
FIGURE 4

Cross-correlation RA(τ) for different lags D with BTc = 10; modulation: BOCs(m, m)

2.1.3 CBOCpilot Signals

Here, the distorted signal is Embedded Image, with the CBOC(6,1,1/11) pilot signal being defined as s(t) = [αp11(t) – βp61(t)]c(t) = αs11(t) – βs61(t), where Embedded Image, Embedded Image, p11(t) and p61(t) are, respectively, the BOCs(1,1) and BOCs(6,1) subcarriers, and c(t) is the code signal. Here, s11(t) = c(t)p11(t) and s61(t) = c(t)p61(t) are, respectively, the BOCs(1,1) and BOCs(6,1) signals (European GNSS Agency, 2023). The receiver correlation scheme is shown in Figure 2, with the local generator given by g(t) = s11(t). As mentioned above, this correlation scheme is different from those of the BPSK and BOCs modulations, as the signal of the local oscillator is not a replica of the nominal transmitted signal, with the purpose of simplifying the receiver architecture. For this case, we have the following:

Embedded Image 17

with Embedded Image (Sousa & Nunes, 2013), where ΛL(|τ| – τ0) = ΛL(τ + τ0) + ΛL(τ – τ0) denotes two triangular pulses of duration 2L centered at τ = ±τ0.

The cross-correlation between signals s61(t) and s11(t) is given by the following (Sousa & Nunes, 2013):

Embedded Image 18

To determine Rϵ,s11(τ), we must account for Equation (1), which is equivalent to ϵ(t) = s(t – d) – s(t) if s(t – d) ≥ s(t) and ϵ(t) = 0 otherwise. Thus, we obtain the following:

Embedded Image 19

Now, using Equation (17) yields Embedded Image. The cross-spectral density corresponding to Embedded Image is as follows:

Embedded Image 20

The cross-correlation at the receiver is given by Equation (7):

Embedded Image 21

with Tc = 1/1.023μs or:

Embedded Image 22

where RA,BOCs(1,1)(τ) is determined by Equation (16) and:

Embedded Image 23

Now, using Equation (18), the cross-spectrum of s11(t) and s61(t) is as follows:

Embedded Image 24

Substituting Equation (24) in Equation (23) yields the following:

Embedded Image 25

2.2 Threat Model TM-A2

In this threat model, the distortion in the subcarriers occurs before modulation by the code signal. This distortion corresponds to the digital distortion 2 described by Pagot et al. (2016), being applicable to signalings in which one or more subcarriers are separated from the code signal, namely, BOCs(m, m) and CBOCpilot signals.

2.2.1 BOCs(m, m) Signals

Let parameter d11 express the subcarrier lead/lag. The distorted subcarrier is Embedded Image, and the error signal is as follows:

Embedded Image 26

Note that for |d11| > Tc/2, the error model in Equation (26) may not be realistic. The modulated signal is s(t) = p11(t)c(t), and its distorted version is Embedded Image, where c(t) is the code signal. The receiver cross-correlation with the local BOCs signal is Embedded Image, which can be written as follows:

Embedded Image 27

where Rs(τ) ≡ 〈s(t)s(t – τ)〉. The cross-correlation Rϵ,11(τ) ≡ 〈ϵ11(t)c(t)s(t – τ)〉 is depicted in Figure 5 and can be expressed by the following difference in triangles:

Embedded Image 28

FIGURE 5
FIGURE 5

Cross-correlation Rϵ,11(τ)

In the frequency domain, Equation (27) can be written as follows:

Embedded Image 29

with Gs(f) given by Equation (15). The cross-spectrum Gϵ,11(f) is determined from Equation (28), leading to the following:

Embedded Image 30

Using Equation (29), we obtain the following:

Embedded Image 31

with D11 ≡ d11/Tc and Tc = 1/1.023μs. Thus, we have the following:

Embedded Image 32

In the absence of EWF distortion, i.e., for D11 = 0, RA(τ) is an even function.

2.2.2 CBOCpilot Signals

Two parameters are required to express the lead/lag of each subcarrier: d11 and d61. These delays are assumed to be independent.

The distorted transmitted signal is Embedded Image, with the distorted subcarriers being Embedded Image for BOCs(1,1) and Embedded Image for BOCs(6,1). The error signals are ϵ11(t), given by Equation (26), and:

Embedded Image 33

Note that for |d61| > Tc/12, the error model in Equation (33) may not be realistic. The receiver cross-correlation with the local BOCs signal is Embedded Image, which can be written as follows:

Embedded Image 34

In the frequency domain, Equation (34) can be written as follows:

Embedded Image 35

with G11(f), G61,11(f), and Gϵ,11(f) given, respectively, by Equations (15), (24), and (30). Moreover, the following can be shown:

Embedded Image 36

leading to:

Embedded Image 37

Accounting for Equation (35) yields the following:

Embedded Image 38

with RA,BOCs(τ) and βC(τ) given, respectively, by Equations (16) and (25). This finally leads to the following result:

Embedded Image 39

In the absence of EWF distortion, i.e., for D61 = D11 = 0, RA(τ) is an even function.

Figure 6 displays the cross-correlation RA(τ) for D11 = 0.15 and D61 = 0.05. The cross-correlation for D11 = D61 = 0 is also shown (non-distorted signal).

FIGURE 6
FIGURE 6

Cross-correlation RA(τ) for TM-A2 with BTc = 8 and different values of D11 and D61; modulation: CBOCpilot

The curves in Figures 3, 4, and 6 show that a code discriminator with early-late spacing ΔEL > 0 will exhibit synchronization errors when the EWF is present because the position of the curves is shifted, in addition to being distorted. This effect will provoke positioning errors in the receiver. Later, we will describe methods that use the curve distortion to detect the presence of EWFs.

3 THREAT MODELS B AND C

TM-B introduces amplitude modulation and model degradations in the analog section of a satellite. More specifically, this model consists of the output from a second-order system when the code-modulated baseband signal is the input (ICAO, 2018; Phelts, 2001). TM-B assumes that the degraded satellite subsystem can be described as a linear system dominated by a pair of complex conjugate poles. The poles are located at σ ± jωd, where σ is the damping factor in nepers/second and fd = ωd/(2π) is the resonant frequency in Hertz. For t ≥ 0, the system response to the step function, with a discontinuity at t = 0, is y(t) = 1 – exp(–σt)[cos(ωdt) + (σ/ωd)sin(ωdt] (D’Azzo et al., 2003). Let ωn and define, respectively, the (undamped) natural frequency and the damping ratio (with ξ < 1), such that σ = ξωn and Embedded Image. The transfer function of the second-order low-pass filter describing the damped linear system is given as follows:

Embedded Image 40

with Ω ≡ 2πf/ωn. Next, we present expressions for the cross-correlation RA(τ) obtained in the scenario corresponding to TM-B:

Embedded Image 41

yielding:

Embedded Image 42

where ρ ≡ 2πB/ωn and μ ≡ ωnTc/(2π).

3.1 BPSK(m) Signals

The power spectrum of BPSK signals is given by Equation (11). Substituting this result in Equation (42) leads to the following:

Embedded Image 43

with:

Embedded Image 44

In Appendix B, it is shown that the cross-correlation RA(τ) can be written as an analytic expression that involves sine and cosine integral functions. Figure 7 illustrates the correlation deformation induced by a TM-B-type distortion. The nominal signal in the figure refers to the cross-correlation in the absence of EWF distortion, which is given in Equation (14) with D = 0.

FIGURE 7
FIGURE 7

Cross-correlation RA(τ) for TM-B with fd = 8 MHz and σ = 8.8 Mnepers/s; modulation: BPSK(1)

In addition to deformation, the EWF also introduces a delay in the cross-correlation, which leads to range errors in the receiver’s DLL. Figure 8 exhibits the range errors in meters for different values of the damped frequency fd = ωd/(2π) and damping factor σ. Note that the largest errors occur when fd and σ are simultaneously small. For simplicity, we assume that the receiver’s DLL tracks the peak of RA(τ).

FIGURE 8
FIGURE 8

Range errors in meters for TM-B with BPSK(1) signals and BTc = 5

3.2 BOCs(m, m) Signals

Using Equation (15), we can derive the following (see also the work by Sousa and Nunes (2013)):

Embedded Image 45

Substituting this result in Equation (42) leads to the following:

Embedded Image 46

where:

Embedded Image 47

Figure 9 illustrates the correlation deformation induced by a TM-B-type distortion. The nominal signal in the figure refers to the cross-correlation in the absence of a TM-B-type distortion, which results from doing D = 0 in Equation (16).

FIGURE 9
FIGURE 9

Cross-correlation RA(τ) for TM-B with fd = 8 MHz and σ = 8.8 Mnepers/s; modulation: BOCs(1,1)

3.3 CBOCpilot Signals

For a CBOCpilot signal affected by TM-B distortion, we can apply Equation (21) with d = 0, which yields the following:

Embedded Image 48

We now incorporate Equations (45) and (24) in Equation (48) and utilize the following definition:

Embedded Image 49

We use the following trigonometric equality:

Embedded Image 50

which allows us to establish the relation between K2 and K0 as follows:

Embedded Image 51

Therefore, we obtain the following:

Embedded Image 52

The cross-correlations RA(τ) of TM-B, defined with the help of functions K0(μ, ξ; x), K1(μ, ξ; x), and K2(μ, ξ; x), can be expressed through analytic formulas using sine and cosine integral functions of complex arguments, as shown in Appendix B for BPSK signals. However, in contrast to the TM-A case, the derivation of cross-correlations for TM-B is generally complicated. An alternative straigh-forward solution resorts to the direct numerical integration of Equations (43), (46), and (52). The IDFT/DFT method, referred to in Section 2, may be used as well.

TM-C introduces both lead/lag and amplitude modulation. Specifically, this model consists of outputs from a second-order system when the modulated baseband signal suffers from lead or lag. This waveform includes a combination of TM-A and TM-B effects (ICAO, 2018). Thus, the main difference between the expressions obtained for TM-A and those to be derived for TM-C is that Equation (6) is now replaced by the following:

Embedded Image 53

where H(f) represents the channel filter (essentially due to the receiver’s front-end filtering effect) and Embedded Image, given by Equation (40), is due to the degradation introduced in the analog section of the satellite.

4 METRICS AND TESTS

Let us consider the scheme shown in Figure 10 for the detection of EWF-distorted GNSS signals constituted by a bank of 2N + 1 correlators, with N > 1 and a correlator spacing equal to Δ. The phase ϵ denotes the code synchronization error. In the work by Pagot et al. (2016), the authors noted that the overall delay interval 2NΔ of the correlator bank need not be equal to twice the chip duration (Tc) for detection purposes. In fact, NΔ ≈ 0.25Tc is generally sufficient because (1) TM-like distortions are more visible around the prompt of the correlation function and (2) correlator outputs situated at too great a distance from the prompt are more subject to multipath.

FIGURE 10
FIGURE 10

Scheme of the correlator bank

The correlator spacing is typically lower-bounded by Δ ≈ 10 ns because the correlator outputs with very small Δ values are strongly correlated (and thus redundant), due to the effect of the front-end filter. Moreover, a time delay of 10 ns between correlators is now reachable, but lower values may be difficult to achieve by analog-to-digital converters.

In Figure 10, r(t) is the (band-pass) received signal, which includes the GNSS signal of power P and AWGN of power spectral density N0/2. The corresponding carrier-to-noise ratio is (C/N0) = P/N0. The baseband signal is Embedded Image, where h(t) is the impulse response of the receiver filter and n(t) is Gaussian low-pass noise with power spectral density Gn(f) = N0|H(f)|2 = N0 if |f| ≤ B and Gn(f) = 0 otherwise. We assume that 2B is the front-end bandwidth and that the receiver’s phase-tracking loop is closed with a sufficient rate to render the effect of the Doppler frequency negligible along the processing chain.

Because 2BT ≫ 1, n(t) is generally well approximated by white noise, that is, E{n(t)n(λ)} ≈ N0δ(t – λ), with E{·} denoting the expectation. The correlator outputs are Gaussian random variables defined as follows:

Embedded Image 54

where the additive noises:

Embedded Image 55

are zero-mean, Gaussian random variables, with the following correlations:

Embedded Image 56

For simulation purposes, the noise vector n = [n–N…n0,…nN]T may be easily generated by using the procedure presented in Appendix C.

The SQM consists of a test to determine whether the signal is affected by EWF distortion. Let us consider the hypotheses H0 (signal is not distorted) and H1 (signal is distorted). Let Embedded Image denote the test variable, which is normalized such that Embedded Image under hypothesis H0 and typically Embedded Image under hypothesis H1. The result of the test is binary: if Embedded Image, it is assumed that the signal is not distorted; otherwise, we assume that the signal is distorted. For a given test, we define the following probabilities:

Embedded Image 57a

Embedded Image 57b

It is not possible to simultaneously minimize Pfa and Pmd. A workable solution to this issue is the Neyman–Pearson criterion, which consists of fixing Pfa at a preselected value and then computing Pmd (Kay, 1998).

For a given metric μ (which is a function of the correlator outputs), the corresponding test variable is defined as follows:

Embedded Image 58

where μmea is the measured value of the metric, which is disturbed by thermal noise and may or may not be affected by signal distortion (anomaly), μnom is the nominal value of the metric (without additive noise or distortion), and λ is a normalization parameter that depends on Pfa. First, λ is adjusted to achieve a desired value of Pfa. Then, Pmd is determined from the values of the correlator measurements. This probability depends on the amount of distortion in the autocorrelation function, the carrier-to-noise ratio, etc.

In work by Pagot et al. (2016), three tests based on ratio metrics were proposed: simple ratio, difference ratio, and sum ratio. These metrics will be analyzed next. The ratio tests (or tests based on ratio metrics) specifically attempt to detect the presence of dead zones (flat correlation peaks) and abnormally sharp or elevated correlation peaks (Phelts, 2001).

4.1 Simple Ratio Metric

Assuming that there are 2N + 1 correlators, it is possible to form 2N simple ratio metrics Embedded Image, i = ±1, …, ± N. The tests are based on random variables expressed by Embedded Image, where Embedded Image. Therefore, Embedded Image.

Under hypothesis H0 (no EWF), the mean and variance of Embedded Image can be approximately calculated using the following second-order Taylor series expansions (Papoulis & Pillai, 2002):

Embedded Image 59

Embedded Image 60

where Embedded Image and Embedded Image are Gaussian random variables with correlation coefficient r = (E{XY} – mXmY)/(σXσY). Considering Equations (54) and (56), we have Embedded Image, Embedded Image, Embedded Image, and r = Rs(iΔ)/Rs(0). Thus, we obtain the following:

Embedded Image 61

When (C/N0)T ≫ 1, we can simplify Equation (61), yielding Embedded Image, which is a convenient result, as it does not depend on the carrier-to-noise ratio. We also have the following:

Embedded Image 62

with Embedded Image. Assuming that Embedded Image is Gaussian under hypothesis H0, we obtain the following probability of false alarm for test Embedded Image:

Embedded Image 63

where erfc(·) denotes the complementary error function:

Embedded Image 64

For the sake of simplicity, let us consider that the 2N parameters Embedded Image are individually adjusted to achieve a common probability of false alarm (Pfa), that is, Embedded Image, i = ±1, …, ± N. The decision criterion consists of declaring an anomaly if at least one of the tests is positive Embedded Image. The corresponding probability of a false alarm is Embedded Image. Using the union bound approximation (Proakis & Salehi, 2008) yields the following:

Embedded Image 65

Taking into account Equations (63) and (65), we obtain the normalization parameter for each pair of correlators:

Embedded Image 66

4.2 Difference Ratio Metric

Considering 2N +1 correlators, i.e., N pairs of early-late correlators, we can form N difference ratio metrics Embedded Image, i = 1,…, N. However, Embedded Image. Thus, the test variables are Embedded Image, with Embedded Image and:

Embedded Image 67

Assuming that Embedded Image is Guassian, the probability of a false alarm for test Embedded Image is given by an expression similar to Equation (63). As in the simple ratio metric, we assume that all tests are characterized by a common probability of false alarm (Pfa), that is, Embedded Image. For each pair of early-late correlators, we obtain a different value of the normalization parameter:

Embedded Image 68

4.3 Sum Ratio Metric

Considering 2N + 1 correlators, we can form N sum ratio metrics Embedded Image, i = 1,…, N. The test variables are Embedded Image, where Embedded Image. In the absence of EWF distortion, the tests have zero mean and variance:

Embedded Image 69

Assuming that Embedded Image is Gaussian, the probabilities of false alarm are expressed by an expression similar to Equation (63), and we obtain a different value of the normalization parameter for each pair of early-late correlators:

Embedded Image 70

5 SIMULATION RESULTS

The ratio tests discussed in the previous section allow for varying performance in EWF detection, depending on the characteristics of the distorted cross-correlation RA(τ). In general, tests based on simple ratio or sum ratio metrics are used to detect deformations in the cross-correlation that result in either increased flatness or sharpness. Tests based on difference ratio metrics excel at detecting asymmetries in the cross-correlation curve but tend to perform poorly when the cross-correlation, although flattened, remains symmetric. The receiver bandwidth also plays an important role in EWF detection. Figure 11 displays the regions of detectability for TM-A1 and BPSK(1) modulation for PFA = 10−4, PMD ≤ 10−3, and a bank of 2N + 1 = 5 correlators with spacing Δ = 0.15Tc. Two different bandwidths were considered: BTc = 2 and 5. The darkest region of each plot indicates the region of detectability. The plots show that the interval of detectability (D ≥ Dmin), defined in terms of the minimum value of parameter D, increases with the value of (C/N0)T in both plots. Additionally, for a given value of (C/N0)T, the interval of detectability grows with the bandwidth.

FIGURE 11
FIGURE 11

Detectability regions for TM-A1 with BPSK(1) and different bandwidths

As noted in Section 2, the presence of a TM-A1 EWF anomaly provokes a deformation and a time shift in the cross-correlation RA(τ) equal to DTc/2, which is responsible for the range error cDTc/2. If we consider that the error for a certain application must not exceed the quantity maximum-allowable error in range (MERR), EWF anomalities with D values less than 2 · MERR/(cTc) may not be detected, but the corresponding range errors will not be hazardous (Rife & Phelts, 2006). Figure 12 displays the hazard regions (darker regions), corresponding to Figure 11, for which the range error exceeds the MERR value of 2 m and the EWF anomaly is not detected by the test.

FIGURE 12
FIGURE 12

Hazard regions for TM-A1 with BPSK(1), MERR = 2 m, and different bandwidths

Figure 13 presents the detectability regions (darker regions) for TM-B with BPSK(1) modulations in the threat space region of 0 < fd < 20 MHz and 0 < σ < 200 Mnepers/s, assuming two signal-to-noise ratios: (C/N0)T = 40 dB and 48 dB. As expected, the detectability region grows with increasing signal-to-noise ratio, starting with small values of fd and σ. To obtain the corresponding hazard regions, for a given MERR value, the plots in Figure 13 could be combined with the range error plot in Figure 8.

FIGURE 13
FIGURE 13

Detectability regions for TM-B with BPSK(1) and different signal-to-noise ratios

6 FINAL REMARKS

In the first part of this work, mathematical expressions were derived for signal cross-correlation functions that are applicable to different threat models and modulations. For TM-A, simple closed-form expressions were derived that can be rapidly computed via sine integral functions with real arguments. For TM-B, only results expressed in terms of integrals could be obtained with generality. When available, specific expressions for the cross-correlation allow for a significant reduction in the computational effort required by SA simulations of different EWF scenarios, with an emphasis on determining the detectability and hazard regions of the threat space. Nevertheless, even when closed-form solutions are not available, numerical integration generally constitutes a rapid computational alternative.

In the second part of this work, a bank of N pairs of early-late (symmetric) correlators combined with a prompt correlator was used with a test based on simple, difference, or sum ratio metrics. Alternatively, we may devise a global test in which the three metrics are used simultaneously, as proposed by Sousa et al. (2022). Although this strategy may introduce some degree of redundancy, it prevents the need to select the optimal test for each type of EWF anomaly, as the nature of the anomalies is usually not known beforehand.

This work permits us to draw the following conclusions:

  1. The expressions obtained herein for calculating the receiver’s cross-correlation are general, providing a useful means for evaluating the performance of current and future EWF detection algorithms via SA simulations. Certain closed-form results assume ideal band-pass filtering, which is considered a good approximation to high-order real filters. The SA results presented herein provide useful insights for designing effective monitoring schemes and examining the sensitivity of receivers to specific threats and errors, without incurring a large computational effort. Note that if the correlations were calculated from signal samples (simple MC simulation), the processing burden would be higher, depending on the sample rate and number of bits per sample.

  2. The EWF detection performance improves with the number of correlator pairs N; however, for N > 4, we did not obtain any significant improvement. However, one must be cautious, as this condition may not hold in all scenarios. Moreover, using equal spacing between correlators may not be the optimal solution.

  3. The regions of detectability tend to increase with the growth of the coherent integration interval T. In performance analyses, the product (C/N0)T is an important parameter to consider. Consequently, for small values of (C/N0), large integration intervals are required. However, even for reasonably large values of T, anomaly detectability may not be guaranteed in significant parts of the threat space. In addition, large values of T imply large delays (latency) in the anomaly detection.

  4. The usefulness of these EWF detection algorithms may be jeopardized in the presence of multipath, as discussed by Nunes and Sousa (2022). Therefore, their applicability to mobile receivers in harsh environments is questionable. Note that fixed reference stations generally permit one to obtain better EWF detection performance than mobile receivers because fixed reference stations are not affected by random dynamic stress and the multipath effect can be made negligible.

  5. Given the limitations of EWF detectors using the SQM approach, the aid of other techniques may be necessary, such as receiver autonomous monitoring techniques, which operate at the navigation solution level.

HOW TO CITE THIS ARTICLE

Nunes, F. D., & Sousa, F. M. G. (2024). Efficient signal quality monitoring of GNSS signals disturbed by evil waveforms. NAVIGATION, 71(4). https://doi.org/10.33012/navi.668

ACKNOWLEDGMENTS

This work was partially supported by FCT-Fundação para a Ciência e Tecnologia, I.P., with project reference UIDB/50008/2020 and DOI identifier https://doi.org/10.54499/UIDB/50008/2020. This research was conducted in part under the European Space Agency project MONINT (Monitoring Solution for New Integrity Applications), with GMVIS, S.A. as the prime contractor and Instituto de Telecomunicações as a subcontractor. The authors express appreciation for the collaboration of A. Negrinho, P. Boto, and P. Fernandes of GMVIS.

A COMPUTATION OF INTEGRAL IS2C

Let us consider the following integral:

Embedded Image 71

For BPSK(m) modulations, we have a = 1, and for BOCs(m,m) modulations, we have a = 1/2. Because sin2 x = (1 – cos(2x)) /2, we can write the following:

Embedded Image 72

We can use the following relation:

Embedded Image 73

to derive:

Embedded Image 74

Embedded Image 75

This result permits us to compute Is2c(a, b; L) in terms of trigonometric and sine integral functions, where the sine integral is defined as follows (Abramowitz & Stegun, 1970):

Embedded Image 76

Because the arguments of the sine integrals in Equation (75) are real, these functions can be computed with very small errors using the polynomial approximations given, for instance, by Abramowitz and Stegun (1970) and Rowe et al. (2015).

B ANALYTIC FORMULA OF RA,BPSK(τ) FOR TM-B

The cross-correlation RA,BPSK(τ) in Equation (43) can be expressed analytically as follows:

Embedded Image 77

where a = πμ, b = 2πμ(τ/Tc), c = 4ξ2 – 2, and:

Embedded Image 78

Embedded Image 79

Embedded Image 80

The quantities Ym are calculated in terms of sine and cosine integral functions of complex argument Si(z) and Ci(z), respectively (see, for instance, Abramowitz and Stegun (1970)):

Embedded Image 81

and:

Embedded Image 82

with:

Embedded Image 83

We also have the following:

Embedded Image 84

and:

Embedded Image 85

with:

Embedded Image 86

The complex numbers ωk, k = 1,…, 4 are the roots of ω4 + cω2 + 1 = 0.

C GENERATION OF ARBITRARY GAUSSIAN VECTORS

Consider the generation of a zero-mean Gaussian noise vector U = GW = [u1,…, uM] with arbitrary covariance matrix C = E{UUT}, using the zero-mean noise vector W = [w1,…, wM] with independent Gaussian components of unity variance. This problem occurs, for instance, in MC simulations. The solution consists of determining G from C = GE{WWT}GT = GGT. Matrix C can be written in terms of its eigenvalues and eigenvectors according to C = VΛVT, where the columns of V are orthonormal eigenvectors and Λ = diag{Λ1, λ2,…, λM} is a diagonal matrix with the eigenvalues of C. Finally, Embedded Image.

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.

REFERENCES

  1. ↵
    1. Abramowitz, M., &
    2. Stegun, I.
    (Eds.). (1970). Handbook of mathematical functions: With formulas, graphs, and mathematical tables. Dover. https://store.doverpublications.com/products/9780486612720?srsltid=AfmBOorS0XzvhT6a0R4KdhNXfgNa7FBJJF5DMJJFapDVfm7Bu0Ht2-0z
  2. ↵
    1. Ali, K.,
    2. Manfredini, E., &
    3. Dovis, F.
    (2014). Vestigial signal defense through signal quality monitoring techniques based on joint use of two metrics. Proc. of the 2013 IEEE/ION Position, Location, and Navigation Symposium (PLANS 2014), Monterey, CA, 1240–1247. https://doi.org/10.1109/PLANS.2014.6851499
  3. ↵
    1. Betz, J. W.
    (2016). Engineering satellite-based navigation and timing: Global navigation satellite systems, signals, and receivers. IEEE Press; Wiley. https://doi.org/10.1002/9781119141167
  4. ↵
    1. D’Azzo, J.,
    2. Houpis, C., &
    3. Sheldon, S.
    (2003). Linear control system analysis and design with MATLAB (6th ed.). CRC Press. https://www.google.com/books/edition/Linear_Control_System_Analysis_and_Desig/R07OBQAAQBAJ
  5. ↵
    1. European GNSS Agency
    . (2023). Galileo open service – signal-in-space interface control document (OS SIS ICD) – issue 2.1. Publications Office of the European Union. https://doi.org/10.2878/39727
  6. ↵
    1. Gómez-Casco, D.,
    2. Crosta, P., &
    3. Spangenberg, M.
    (2022). Validation of evil waveforms in a GNSS simulator for GPS and Galileo signals. Proc. of the 10th Workshop on Satellite Navigation Technology (NAVITEC 2022), Noordwijk, Netherlands, 1–7. https://doi.org/10.1109/NAVITEC53682.2022.9847541
  7. ↵
    1. ICAO
    . (2018). Annex 10 – Aeronautical Telecommunications – Volume 1 – Radio Navigation Aids (7th ed.). International Civil Aviation Organization. https://store.icao.int/en/annex-10-aeronautical-telecommunications-volume-i-radio-navigational-aids
  8. ↵
    1. Irsigler, M.
    (2008). Multipath propagation, mitigation and monitoring in the light of Galileo and the modernized GPS [Doctoral dissertation, Bundeswehr University, Munich, Germany]. https://athene-forschung.unibw.de/doc/86276/86276.pdf
  9. ↵
    1. Kay, S.
    (1998). Fundamentals of statistical signal processing: Detection theory. Vol. II. Prentice-Hall. https://www.pearson.com/en-us/subject-catalog/p/fundamentals-of-statistical-signal-processing-detection-theory-volume-2/P200000009354/9780135041352
  10. ↵
    1. Negrinho, A.,
    2. Fernandes, P.,
    3. Boto, P.,
    4. Nunes, F., &
    5. Sousa, F.
    (2021). Evil waveforms detection solutions for autonomous navigation integrity. Proc. of the 34th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2021), St. Louis, MO, 4097–4115. https://doi.org/10.33012/2021.18070
  11. ↵
    1. Nunes, F., &
    2. Sousa, F.
    (2022). Performance of GNSS evil waveform detectors in the presence of multipath. Proc. of the 10th Workshop on Satellite Navigation Technology (NAVITEC 2022), Noordwijk, Netherlands, 1–10. https://doi.org/10.1109/NAVITEC53682.2022.9847550
  12. ↵
    1. Oppenheim, A., &
    2. Schafer, R.
    (1975). Digital signal processing. Prentice-Hall. https://books.google.com/books/about/Digital_Signal_Processing.html?id=vSzuVLBbp6cC
  13. ↵
    1. Pagot, J.
    (2016). Modeling and monitoring of new GNSS signal distortions in the context of civil aviation [Doctoral dissertation, MITT, Toulouse, France]. https://theses.hal.science/tel-01528481
  14. ↵
    1. Pagot, J.,
    2. Julien, O.,
    3. Thevenon, P.,
    4. Julien, O.,
    5. Fernandez, F. A., &
    6. Cabantous, M.
    (2018). Signal quality monitoring for new GNSS signals. NAVIGATION, 65(1), 83–97. https://doi.org/10.1002/navi.218
  15. ↵
    1. Pagot, J.,
    2. Thevenon, P.,
    3. Julien, O.,
    4. Fernandez, F. A., &
    5. Maillard, D.
    (2016). Threat models design for new GNSS signals. Proc. of the 2016 International Technical Meeting of the Institute of Navigation, Monterey, CA, 970–982. https://doi.org/10.33012/2016.13476
  16. ↵
    1. Papoulis, A., &
    2. Pillai, S. U.
    (2002). Probability, random variables, and stochastic processes (4th ed.). McGraw-Hill. https://www.google.com/books/edition/_/cUmiDAEACAAJ
  17. ↵
    1. Phelts, R.
    (2001). Multicorrelator techniques for robust mitigation of threats to GPS signal quality [Doctoral dissertation, Stanford University]. https://web.stanford.edu/group/scpnt/gpslab/pubs/theses/EricPheltsThesis01.pdf
  18. ↵
    1. Phelts, R., &
    2. Akos, D.
    (2006). Effects of signal deformations on modernized GNSS signals. Journal of Global Positioning Systems, 5(2), 2–10. https://web.stanford.edu/group/scpnt/gpslab/pubs/papers/Phelts_JournalGPS_2006.pdf
  19. ↵
    1. Pirsiavash, A.,
    2. Broumandan, A.,
    3. Lachapelle, G., &
    4. O’Keefe, K.
    (2017). Detection and classification of GNSS structural interference based on monitoring the quality of signal at the tracking level. Proc. of the 6th ESA International Colloquium on Scientific and Fundamental Aspects of the Galileo, Valencia, Spain, 1–10. https://www.researchgate.net/publication/321263970_Detection_and_Classification_of_GNSS_Structural_Interference_Based_on_Monitoring_the_Quality_of_Signal_at_the_Tracking_Level
  20. ↵
    1. Proakis, J., &
    2. Salehi, M.
    (2008). Digital communications (5th edition). McGraw-Hill. https://www.google.com/books/edition/Digital_Communications/ABSmAQAACAAJ
  21. ↵
    1. Raimondi, M.,
    2. Sénant, E.,
    3. Fernet, C.,
    4. Pons, R., &
    5. AlBitar, H.
    (2012). Generating evil waveforms on Galileo signals using NAVYS. Proc. of the 6th ESA Workshop on Satellite Navigation Technologies (Navitec 2012) & European Workshop on GNSS Signals and Signal Processing, Noordwijk, Netherlands, 1–8. https://doi.org/10.1109/NAVITEC.2012.6423071
  22. ↵
    1. Rife, J., &
    2. Phelts, E.
    (2006). Formulation of a time-varying maximum allowable error for ground-based augmentation systems. Proc. of the 2006 National Technical Meeting of the Institute of Navigation, Monterey, CA, 441–453. https://www.ion.org/publications/abstract.cfm?articleID=6550
  23. ↵
    1. Rowe, B. T.,
    2. Jarvis, M.,
    3. Mandelbaum, R.,
    4. Bernstein, G. M.,
    5. Bosch, J.,
    6. Simet, M.,
    7. Meyers, J. E.,
    8. Kacprzak, T.,
    9. Nakajima, R.,
    10. Zuntz, J.,
    11. Miyatake, H.,
    12. Dietrich, J. P.,
    13. Armstrong, R.,
    14. Melchior, P., &
    15. Gill, M. S. S.
    (2015). GalSim: The modular galaxy image simulation toolkit. Astronomy and Computing, 10, 121–150. https://doi.org/10.1016/j.ascom.2015.02.002
  24. ↵
    1. Selmi, I.,
    2. Thevenon, P.,
    3. Macabiau, C.,
    4. Olivier, J., &
    5. Mabilleau, M.
    (2020). Signal quality monitoring algorithm applied to Galileo signals for large evil waveform threat space. Proc. of the 2020 International Technical Meeting of the Institute of Navigation, San Diego, CA, 352–365. https://doi.org/10.33012/2020.17149
  25. ↵
    1. Shallberg, K.,
    2. Ericson, S.,
    3. Phelts, E.,
    4. Walter, T.,
    5. Kovach, K., &
    6. Altshuler, E.
    (2017). Catalog and decription of GPS and WAAS L1 C/A signal deformation events. Proc. of the 2017 International Technical Meeting of the Institute of Navigation, Monterey, CA, 508–520. https://doi.org/10.33012/2017.14877
  26. ↵
    1. Sousa, F., &
    2. Nunes, F.
    (2013). New expressions for the autocorrelation function of BOC GNSS signals. NAVIGATION, 60(1), 1–9. https://doi.org/10.1002/navi.30
  27. ↵
    1. Sousa, F.,
    2. Nunes, F.,
    3. Negrinho, A.,
    4. Fernandes, P., &
    5. Boto, P.
    (2022). Signal quality monitoring aspects in GNSS signals affected by evil waveforms. Proc. of the 10th Workshop on Satellite Technology (NAVITEC 2022), Noordwijk, Netherlands, 1–14. https://doi.org/10.1109/NAVITEC53682.2022.9847553
  28. ↵
    1. Thevenon, P.,
    2. Julien, O.,
    3. Tessier, Q.,
    4. Maillard, D.,
    5. Cabantous, M.,
    6. Amarillo-Fernández, F., &
    7. Oliveira-Salgueiro, F.
    (2014). Detection performances of evil waveform monitors for the GPS L5 signal. Proc. of the 27th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2014), Tampa, FL, 3331–3342. https://www.ion.org/publications/abstract.cfm?articleID=12453
  29. ↵
    1. Tranter, W. H.,
    2. Shanmugan, K. S.,
    3. Rappaport, T. S., &
    4. Kosba, K. L.
    (2004). Principles of communication systems simulation with wireless applications. Prentice-Hall. https://dl.acm.org/doi/book/10.5555/1407789
Back to top

In this issue

Print
Download PDF
Article Alerts
Email Article
Citation Tools
Efficient Signal Quality Monitoring of GNSS Signals Disturbed by Evil Waveforms
Fernando D. Nunes, Fernando M. G. Sousa
NAVIGATION: Journal of the Institute of Navigation Dec 2024, 71 (4) navi.668; DOI: 10.33012/navi.668
Twitter logo Facebook logo Mendeley logo
Bookmark this article

Jump to section

Keywords