GNSS interference detection: Methodology utilising ADS-B NACp indicator and GPS almanac data

Abstract

The threat of GNSS interference poses a great danger to many critical infrastructure systems including air navigation. With a focus on mitigating this threat, this paper proposes a methodology for detecting GNSS interference. The methodology utilises the quality indicator NACp transmitted in ADS-B messages and GPS almanac data for interference detection. The NACp indicator enables estimation of the position error derived from GPS, which is compared with the HDOP value of the GPS satellite constellation. Based on this comparison, the developed detection algorithm determines whether the aircraft is affected by jamming. The detection methodology is evaluated on datasets obtained during deliberate experiments with GPS jamming. The proposed methodology provides a way to detect GNSS interference, facilitating mitigation of its impact on air traffic operation.

Nomenclature

${E_k}$

eccentric anomaly

$Q$

covariance matrix

${R_c}$

containment radius

ACC

accuracy

ADS-B

Automatic Dependent Surveillance-Broadcast

ADS-B 1090 ES

ADS-B 1090 MHz extended squitter

API

application programming interface

ATM

air traffic management

CTU

Czech Technical University in Prague

ECEF

Earth-Centred, Earth-Fixed

EPU

estimated position uncertainty

FN

false negatives

FP

false positives

FPR

false positives rate

GNSS

Global Navigation Satellite System

GPS

Global Positioning System

HDOP

horizontal dilution of precision

HFOM

Horizontal Figure of Merit

HIL

Horizontal Integrity Limit

HPL

Horizontal Protection Level

ICAO

International Civil Aviation Organization

MISC

misclassification rate

MOPS

Minimum Operational Performance Standards

NANU

Notice Advisory to Navstar Users

NACp

Navigation Accuracy Category for Position

NIC

Navigation Integrity Category

PPV

Precision

RFI

Radio Frequency Interference

RTCA

Radio Technical Commission for Aeronautics

SA

selective availability

SBAS

Satellite Based Augmentation Systems

TN

true negatives

TP

true positives

TPR

true positives rate

UERE

User Equivalent Range Error

UTC

Universal Time Coordinated

WGS

World Geodetic System

Greek symbol

$\sigma $

error of the measured pseudo-range

1.0 Introduction

The resilience of communication, navigation and surveillance systems is essential to ensure the safety of air traffic operations. As aviation relies on Global Navigation Satellite Systems (GNSS) as the primary navigation source, it is important to detect GNSS radio frequency interference (RFI) and to smoothly transition to complementary navigation capabilities provided by ground-based systems such as distance measuring equipment [1, 2]. In addition to GNSS RFI detection, the technology providing resilience to interference is also an important component in mitigating the negative effects of GNSS RFI. Resilient receivers, active antenna arrays capable of adapting their radiation patterns, aviation dual-frequency, multi-constellation receivers and other resilient technologies are being developed by multiple research teams [3–5].

Several working papers highlighting the increasing concerns about GNSS vulnerabilities, including RFI have been presented in global aviation meetings of the International Civil Aviation Organization (ICAO) in recent years [2, 6, 7]. These working papers call for implementing countermeasures to improve GNSS resilience in both the short and long term. These countermeasures shall preserve the safety, capacity and environmental benefits of air traffic management (ATM) and its key GNSS-based tools like the performance-based navigation concept and the Automatic Dependent Surveillance-Broadcast (ADS-B) surveillance technology.

The EUROCONTROL Think Paper #9 [8] brought an assessment of the operational impact of the threat of RFI to satellite navigation in aviation, by presenting data derived from the EUROCONTROL Voluntary ATM Incident Reporting System. The document shows that the impact of GNSS RFI is particularly severe in conflict zones. However, RFI is not limited to these areas, since up to 38.5% of European en-route traffic areas are affected by regular GNSS RFI. This increases the workload of pilots and air traffic controllers thus causing a significant safety risk. Additionally, GNSS interference monitoring conducted on European airports shows regular occurrence of GNSS RFI in their vicinity [9, 10].

The first step in mitigating the impact of GNSS RFI is its detection which enables taking appropriate actions to ensure a smooth continuation of air traffic. For this purpose, the EUROCONTROL Handbook No. 6 [1] has been published in 2021. The Handbook includes samples of contingency scenarios, including the case of GPS loss due to intentional interference in the form of jamming. For timely activations of these scenarios, GNSS jamming detection is a key component.

One of the ways to detect GNSS jamming is by using Automatic Dependent Surveillance-Broadcast 1090 MHz Extended Squitter (ADS-B 1090 ES) surveillance technology. Several studies have been conducted using quality indicators from ADS-B messages to detect GNSS interference [11–20]. The most common ADS-B 1090 ES parameters used for the detection of GNSS RFI are the following:

• • Navigation Accuracy Category for Position (NACp) – provides an indication of the position accuracy based on the Estimated Position Uncertainty (EPU) parameter. EPU is defined as the radius of a circle, centred at the calculated position on an aircraft. This circle contains the actual aircraft position at least 95% of the time [21, 22]. When the position information is generated by a GPS or GNSS system, EPU is also commonly referred to as a Horizontal Figure of Merit (HFOM).

• • Navigation Integrity Category (NIC) – provides information about the value of an integrity containment radius ( ${R_c}$ ). This value is derived from the Horizontal Protection Level (HPL) or the Horizontal Integrity Limit (HIL) [21].

An example of a detection methodology using only the NIC indicator has been proposed by Standford University [14]. Several follow-up papers focused on determining the size and shape of the area affected by RFI as well as estimating the source of the interference were produced by the same research group [13, 15, 16]. These algorithms typically require ADS-B data collected over several hours with the presence of long-term jamming. Another detection method using the NIC indicator was put forward by EUROCONTROL [11, 19], this method aims at identifying the source of the interference. Alternatively, the NACp indicator can also be used for interference detection as shown on ADS-B data from the OpenSky Network [23]. Various approaches to GNSS RFI detection using pattern recognition and cluster analysis [17], leveraging the correlation between NACp and NIC indicators [18], or constructing discretestatistical mathematical model [20] have been developed at the Czech Technical University in Prague (CTU). Some research also focused on detecting GPS spoofing from ADS-B [12], however spoofing detection is out of the scope of the presented methodology.

The objective of the methodology presented is to provide an additional or complementary means of determining whether a GNSS signal received by an aircraft is affected by intentional jamming. The methodology utilises data transmitted by aircraft via the ADS-B 1090 ES technology which is further specified in ICAO Annex 10 Volume IV, Surveillance Radar and Collision Avoidance Systems [24], and the Minimum Operational Performance Standards for ADS-B 1090 ES [21, 22]. ADS-B data is freely available to all entities.

The proposed methodology only considers ADS-B version 2 standard. Version 2 is obligatory in the air segment of the European Region according to European Commission Regulation (EU) 1207/2011 [25] and its subsequent amendments. In June 2023, the equipage of ADS-B version 2 is 95.2% for aircraft and 97.3% for flights in the European area operated by EUROCONTROL [26]. The methodology does not cover the 1090 ES ADS-B technologies of ICAO versions 0 and 1, nor the functionalities and parameters available in future ICAO version 3.

In general, there are two cases of GNSS jamming interference affecting an aircraft and therefore ADS-B. In the first case, a complete loss of positional information occurs, while in the second case, the positional information is only degraded. The majority of the aforementioned research on GNSS RFI detection from ADS-B focuses only on the complete loss of positional information. Additionally, this loss must impact multiple aircraft to ensure detection of the jamming. Conversely, the proposed methodology focuses on detecting short-term interference as well as position information degradation, which is more challenging to detect compared to long-term jamming which impacts multiple aircraft.

The methodology shows that GNSS interference can be detected by examining real-time changes in the NACp indicator and comparing them to a pessimistic estimate of the position inaccuracy due to the position of GPS satellites. Thus, the method enables real-time interference detection capabilities through ADS-B. This enables air navigation service providers to be informed about the presence of jamming in a timely manner which subsequently allows them to prepare airspace users for potential GNSS degradation or unavailability, increasing aviation safety in the process.

One of the big advantages of this method is the high accessibility of both the ADS-B and GPS almanac data which greatly simplifies its implementation. In comparison to the previously published GNSS RFI detection methods, a minimum of two messages, one affected by jamming and one not affected by jamming, are necessary for jamming detection. In other words, it is not necessary to have a large amount of data from a GNSS jamming event. The methodology can be used independently or in combination with existing detection methods.

The methodology describes procedures for determining the current horizontal position error from satellites and explains the comparison with the NACp value, which leads to the determination of the impact of jamming on an aircraft. In the next sections, the detection results are evaluated using data from deliberate GNSS jamming experiments. The discussion and conclusion sections present the findings and highlight the possible use of the research results.

2.0 Methodology

The methodology is based on determining the maximum possible actual change in the reported position error through the NACP indicator. It estimates whether a decrease in this indicator could occur due to standard effects and conditions, or whether this decrease can no longer be explained by natural causes. The most pessimistic natural estimate of the change in the NACp indicator is derived from the horizontal dilution of precision (HDOP) parameter. The algorithm implementing the methodology is written in Python programming language.

2.1 Calculation of HDOP value from GPS almanacs

The value of HDOP can be calculated by knowing the position of the receiver (aircraft) and by knowing the position of the satellites used in the calculation of the receiver’s position. The methodology assumes that the aircraft GNSS receiver uses only GPS constellation. The positions of the GPS satellites can be calculated at any given moment using the GPS Almanac [27].

However, which satellites are actually used in the position calculation cannot be so clearly determined. Factors like the individual elevation mask for a particular receiver determine what satellites with low elevation angles are used. Therefore, a certain uncertainty has to be allowed in the HDOP estimation. To resolve this, the methodology computes the most pessimistic possible configuration of GPS satellites, resulting in the most pessimistic estimate of HDOP. For low-flying aircraft, satellite visibility can be obstructed by mountainous terrain. To take this factor into account, the methodology implements a digital elevation model of the terrain when calculating the HDOP parameter.

The input to the detection algorithm is ADS-B position message data, which includes the message reception timestamp, latitude, longitude and altitude of the aircraft. As an additional input, altitude in metres may be included to account for terrain obstacles. This is recommended for low-flying aircraft and aircraft in mountainous regions to exclude satellites that are not visible to the aircraft due to terrain obstacles.

The reception timestamp of the ADS-B message plays an important role in the selection of almanac data. The GPS almanac is available from the CelesTrak website [27]. The Yuma almanac format was selected for the algorithm. The Yuma Almanac is published regularly seven times a week, always at the same time. Another important factor is the time difference between the UTC time of the received ADS-B message, where the leap seconds are corrected, and the GPS time. This time difference must be taken into account in order to select the correct Yuma almanac. The algorithm selects the almanac data closest to the date and time of the received ADS-B message. Once the correct Yuma almanac is selected, the values necessary for real-time position error calculation are derived from it. The values obtained from the GPS almanac data are summarised in Table 1.

Table 1. Parameters for each GPS satellite obtained from the Yuma almanac [28]

Knowing the desired time at which the HDOP value shall be determined, the position of all GPS satellites is calculated using the equations in Table 2. The resulting coordinates of the GPS satellites are in the Earth-centred, Earth-fixed (ECEF) coordinate system.

Table 2. Equations used to calculate the position of GPS satellites [28]

The value of eccentric anomaly ${E_k}$ used in one of the equations in Table 2 is found iteratively as a root of the non-linear equation using the scipy library function fsolve [29] based on the Powell’s hybrid method [30]. The constants necessary for calculating the position of GPS satellites and their pseudoranges are as follows:

• • Speed of light [m/s], $c = {\rm{\;\;}}299792458{\rm{\;\;}}\frac{m}{s}$

• • Earth rotation angular velocity, ${\omega _E} = 7.2921151467 \times {10^{ - 5}}\frac{{rad}}{s}$

• • Earth geocentric gravitational constant, $\mu = 398600.5 \times {10^8}{\rm{\;\;}}\frac{{{m^3}}}{{{s^2}}}$

• • Earth inclination angle, $i = 0.3 \times \pi {\rm{\;\;}}rad$

After determining the position of all GPS satellites, the next step is to determine the visible satellites. Starting with excluding satellites with a health status different to “000”, meaning any malfunctioning satellites are disregarded. Second, satellites below the elevation mask angle, which is set to 5° by default [31], are excluded. Finally, for aircraft flying under a certain altitude, terrain obstacles between the aircraft and the satellite are checked. The altitude checking threshold is set to 9000 m, taking into account the highest mountain peak on Earth.

The terrain check is conducted using the Open-Elevation API [32] as a source of terrain elevation data. The distance between the satellite and the aircraft is divided into a number of points. For each point, the elevation is obtained from the Open-Elevation API [32] and it is checked whether the line connecting the aircraft and the satellite is intersected by the terrain. If there is an intersection, the satellite is excluded from the HDOP calculation. An example of aircraft flying in a fjord is seen in Fig. 1.

Figure 1. An example of terrain obstacle between aircraft and satellite with an elevation angle of 16.7°. In this case, the orange aircraft-satellite line of sight intersects the terrain in blue, hence the satellite is excluded from the HDOP calculation.

In the next step, when all visible satellites are determined, the HDOP value at the aircraft location is determined. The HDOP value can be calculated from the covariance matrix $Q$ . An example of a calculation for four satellites is shown in Equations (1), (2), and (3).

(1) \begin{align} G = \left( \begin{array}{l@{\quad}l@{\quad}l@{\quad}l} {{\rm{cos}}\,{e_1}\,{\rm{cos}}\,{a_1}} & {}{{\rm{cos}}\,{e_1}\,{\rm{sin}}\,{a_1}} & {}{{\rm{sin}}\,{e_1}} & {}{ - 1} \\[3pt] {{\rm{cos}}\,{e_2}\,{\rm{cos}}\,{a_2}} & {}{{\rm{cos}}\,{e_2}\,{\rm{sin}}\,{a_2}} & {}{{\rm{sin}}\,{e_2}} {}& { - 1}\\[3pt] {{\rm{cos}}\,{e_3}\,{\rm{cos}}\,{a_3}} & {}{{\rm{cos}}\,{e_3}\,{\rm{sin}}\,{a_3}} & {}{{\rm{sin}}\,{e_3}} & {}{ - 1}\\[3pt] {{\rm{cos}}\,{e_4}\,{\rm{cos}}\,{a_4}} & {}{{\rm{cos}}\,{e_4}\,{\rm{sin}}\,{a_4}} & {}{{\rm{sin}}\,{e_4}} & {}{ - 1} \end{array} \right) = \left( \begin{array}{l@{\quad}l@{\quad}l@{\quad}l}{ - \frac{{x - {x_1}}}{{{r_1}}}} & {}{ - \frac{{y - {y_1}}}{{{r_1}}}} & {}{ - \frac{{z - {z_1}}}{{{r_1}}}} & {}{ - 1}\\[3pt] { - \frac{{x - {x_2}}}{{{r_2}}}} & {} { - \frac{{y - {y_2}}}{{{r_2}}}} & {}{ - \frac{{z - {z_2}}}{{{r_2}}}} & {}{ - 1}\\[3pt] { - \frac{{x - {x_3}}}{{{r_3}}}} & {}{ - \frac{{y - {y_3}}}{{{r_3}}}} & {}{ - \frac{{z - {z_3}}}{{{r_3}}}} & {}{ - 1}\\[3pt] { - \frac{{x - {x_4}}}{{{r_4}}}} & {}{ - \frac{{y - {y_4}}}{{{r_4}}}} & {}{ - \frac{{z - {z_4}}}{{{r_4}}}} & {}{ - 1} \end{array} \right)\end{align}

(2) \begin{align} Q = {\left( {{G^T}G} \right)^{ - 1}} \times {\sigma ^2} = \left( \begin{array}{l@{\quad}l@{\quad}l@{\quad}l} {{d_{11}}} & {}{{d_{12}}} & {}{{d_{13}}} & {}{{d_{14}}}\\[3pt] {{d_{21}}} & {}{{d_{22}}} & {}{{d_{23}}} & {}{{d_{24}}}\\[3pt] {{d_{31}}} & {}{{d_{32}}} & {}{{d_{33}}} & {}{{d_{34}}}\\[3pt] {{d_{41}}} & {}{{d_{41}}} & {}{{d_{43}}} & {}{{d_{44}}} \end{array} \right) \times {\sigma ^2} = \left( \begin{array}{l@{\quad}l@{\quad}l@{\quad}l}{\sigma _x^2} & {}{{\sigma _{xy}}} & {}{{\sigma _{xz}}} & {}{{\sigma _{xt}}}\\[3pt] {{\sigma _{yx}}} & {}{\sigma _y^2} & {}{{\sigma _{yz}}} & {}{{\sigma _{yt}}}\\[3pt] {{\sigma _{zx}}} & {}{{\sigma _{zy}}} & {}{\sigma _z^2} & {}{{\sigma _{zt}}}\\[3pt] {{\sigma _t}} & {}{{\sigma _t}} & {}{{\sigma _t}} & {}{\sigma _t^2} \end{array} \right) \end{align}

(3) \begin{align} HDOP = {\rm{\;\;}}\sqrt {{d_{11}} + {d_{22}}} = {\rm{\;\;}}\frac{1}{\sigma }\sqrt {\sigma _x^2 + \sigma _y^2.} \end{align}

Where, $\sigma $ is the measured pseudorange error (UERE – User Equivalent Range Error), ${e_n}$ elevation angle of the satellite, ${a_n}$ azimuth to the satellite, $x,y,z$ aircraft coordinates in ECEF, ${x_n},{y_n},{z_n}$ coordinates of the satellite in ECEF, and ${r_n}$ – range between the aircraft and the satellite.

2.2 HFOM value transmitted by aircraft

The ADS-B quality indicator NACp is determined from the HFOM parameter transmitted together with the position data from the aircraft’s GNSS receiver. The relation between HFOM (also referred to as EPU) and the NACp indicator is shown in Table 3.

Table 3. Relation between HFOM (also referred to as EPU) and the NACp indicator [21]

HFOM represents the horizontal navigation accuracy at the 95% probability level, assuming a fail-safe condition at the time of use. According to the Radio Technical Commission for Aeronautics’s document DO-229 [33], HFOM is the radius of a circle in the horizontal plane, specifically the local plane tangent to WGS-84 ellipsoid, centred at the actual aircraft position. The way HFOM is calculated may differ based on the used model of a GNSS receiver. In general, the determination of the HFOM value is dependent on two variables, HDOP and $\sigma $ , as shown in Equation (4) [33].

(4) \begin{align} HFOM = 2 \times \sqrt {{\sigma _x}^2 + {\sigma _y}^2} = 2 \times \frac{1}{\sigma }\sqrt {{\sigma _x}^2 + {\sigma _y}^2} \times \sigma = 2 \times \sqrt {{d_{11}} + {d_{22}}} \times \sigma = 2 \times HDOP \times \sigma \end{align}

2.3 Differences in GNSS receivers in aviation

In Equation (4), ${\sigma _x}^2$ and ${\sigma _y}^2$ represent elements on the diagonal of the covariance matrix $Q$ of the position error vector ${(x,y)^T}$ . $\sigma $ in Equation (4) is a theoretical error of the measured pseudo-range at the probability level of $1\sigma = 68{\rm{\% }}$ , assuming the error is the same for all satellites. The sigma parameter is determined differently based on the GNSS receiver type. In general, there are three types of aircraft receivers in aviation, differing in selective availability (SA) awareness and GPS augmentation:

1. 1. GPS receivers without augmentation, not aware of SA off

2. 2. GPS receivers without augmentation, aware SA off

3. 3. GPS receivers augmented with satellite-based augmentation systems (SBAS)

The GNSS RFI detection algorithm differs for GNSS receivers without SBAS augmentation (types 1 and 2) and for SBAS-augmented GNSS receivers (type 3). As far as type 1 receivers are concerned, there may still be very old installations where the GPS receiver works with the existence of SA. The GNSS receiver type is not possible to 100% determine unless provided by the aircraft. However, it can be estimated with some degree of certainty. The type of GNSS receiver is derived from the value of the NACp indicator transmitted by a given aircraft. The transmitted NACp indicator also determines how some parameters in the algorithm are calculated.

Given the increase in accuracy provided by SBAS corrections, it can be assumed that once the aircraft transmits the value of NACp greater than 9, it is considered to be using an SBAS-augmented receiver type 3. In case all the NACp values transmitted by a given aircraft are equal or lower than 9, i.e. NACp $ \le $ 9, it is assumed that the aircraft is equipped with a GNSS type 1 or type 2 receiver. In these types of receivers, the $\sigma $ value should not change rapidly in the short term since the variability of $\sigma $ results from the application of ionospheric corrections according to the Klobuchar model. Without SBAS, the variability is not necessary to consider.

The Klobuchar model does not represent a real-time ionosphere stage. Therefore parameters of the model can change at intervals of one to six days, and all satellites use the same values of the Klobuchar model parameters. It is assumed that GNSS receivers type 1 and 2 should transmit NACp equal to 8 or 9 under standard conditions, meaning without the presence of RFI. On the assumption that the estimate of the pseudorange error has a time-invariant standard deviation $\sigma $ , the maximum $\sigma $ value can be refined using the changing value of the HDOP over a longer period. Therefore, the GNSS RFI detection algorithm uses data from longer time intervals for non-SBAS receivers and current values for the SBAS-augmented type 3 receivers.

2.4 Parameters used in the interference detection algorithm

The detection algorithm evaluates whether the position inaccuracy represented by the NACp indicator could have been caused by standard effects and conditions, like those mentioned in Section 2.1 or whether it is caused by GNSS RFI. In addition to this, the algorithm also determines when the aircraft is no longer affected by the interference. Several intermediate calculations are made within the algorithm from the ADS-B and GPS almanac data. In each loop of the algorithm, the following inputs must be obtained from ADS-B messages:

• • Timestamp

• • Aircraft latitude

• • Aircraft longitude

• • Aircraft altitude

• • NACp indicator

To run the algorithm, it is necessary to receive at least two ADS-B messages. When the first message from an aircraft is received, it is presumed not to be affected by any interference. A new loop of the algorithm is executed after receiving every new ADS-B message. For future reference, the previous loop and its message will be referred to as $t - 1$ and the current loop and message as $t$ . In each loop, the following parameters can be calculated:

• • $HFO{M_{max}}$

• • $HDOP$

• • ${\sigma _{max}}$ , respectively ${\hat \sigma _{max}}$

• • $ HDO{P_{pessimistic}} $

• • $HFO{M_{pessimistic}}\!\left( t \right)$ , respectively $\widehat{HFOM}_{pessimistic}\!\left( t \right)$

• • $ NAC{p_{min}} $

• • $NAC{p_{ref}}$ , respectively $\widehat{NACp}_{ref}$

The following paragraphs describe the calculation of each of the parameters computed within the detection algorithm. These parameters are necessary for the interference detection process.

1. 1. The $HFO{M_{max}}$ value, as the maximum possible position error in meters, is determined from the NACp indicator value using conversion depicted in Table 3. For example, if the received $NACp = 8$ , then the $HFO{M_{max}} = 92.6$ meters.

2. 2. HDOP is calculated using the aircraft position data and the positions of all satellites, as described in Section 2.1. The theoretical best value of HDOP is estimated by means of Equation (3).

   In the case of real-time evaluation, it is necessary to provide satellite status data to exclude unhealthy satellites from the calculation. As an alternative to the satellite health parameter described in Section 2.1, it is possible to use SBAS Message Type 6 [34]. Another way is to use the NANU messages published jointly by the United States Coast Guard and the GPS Operations Center [35]. In this paper, the health of satellites was determined by NANU messages.

3. 3. Having the $HFO{M_{max}}$ and the HDOP values, the most pessimistic sigma value denoted as ${\sigma _{max}}$ is calculated using Equation (5).

   (5) \begin{align} {\sigma _{max}} = \frac{{HFO{M_{max}}}}{{2 \times HDOP}} \end{align}

4. 4. Since the aircraft does not use all theoretically trackable satellites, for example, due to the application of the elevation mask, the HDOP value computed in point 2 must be adjusted. The most pessimistic denoted as $HDO{P_{pessimistic}}\left( t \right)$ be determined as $HDO{P_{pessimistic}}\left( t \right) = max\left\{ {HDOP,1.25} \right\}$ . This is possible since the HDOP value is with 95% probability greater than 1.25 for the GPS constellation [36]. The probability distribution function of HDOP is displayed in Fig. 2.

   

   Figure 2. The cumulative distribution function of HDOP for GPS [36].

5. 5. The largest possible position error caused by standard effects and conditions is denoted as $\widehat{HFOM}_{pessimistic}\!\left( t \right)$ for GNSS receivers type 1 and 3, and as $HFO{M_{pessimistic}}\!\left( t \right)$ for receivers type 3. This value is crucial for decision-making whether the aircraft is affected by jamming. Its calculation differs based on the type of GNSS receiver.

   1. 5A. For receivers type 1 and 2, the $\widehat{HFOM}_{pessimistic}\!\left( t \right)$ value is determined using Equation (6):

      (6) \begin{align} \widehat{HFOM}_{pessimistic}\!\left( t \right) = 2 \times HDO{P_{pessimistic}}\left( t \right) \times {\hat \sigma _{max}}\left( {t - {t_n}} \right) \end{align}

      Where ${\hat \sigma _{max}}\left( {t - {t_n}} \right)$ value represents the minimal value from the time interval when the aircraft was evaluated as non-jammed. The ${t_n}$ refers to the last loop when the aircraft was not affected by jamming. For receivers type 1 and 2, the ${\hat \sigma _{max}}$ is defined as the minimum value of ${\sigma _{max}}$ on the time interval $\langle {t_1},{t_2}\rangle $ . This time interval is defined by the first message received from the aircraft ${t_1}$ and the first message evaluated as jammed ${t_2}$ . In the case of previous jamming impact, the $\langle {t_1},{t_2}\rangle $ time interval is the period during which the aircraft was evaluated as not being jammed between two jamming occurrences.

   2. 5B. For type 3 GNSS reveivers, the $HDO{P_{max}}\left( t \right)$ value is determined using Equation (7):

      (7) \begin{align} HFO{M_{pessimistic}}\!\left( t \right) = 2 \times HDO{P_{pessimistic}}\left( t \right) \times {\sigma _{max}}\left( {t - {t_n}} \right) \end{align}

      Where ${\sigma _{max}}\left( {t - {t_n}} \right)$ represents the last estimated ${\sigma _{max}}$ value when the arcraft was evaluated as not jammed. For receivers type 3, no estimation based on time interval is applied.

6. 6. Next, the calculated $\widehat{HFOM}_{pessimistic}\!\left( t \right)$ or $HFO{M_{pessimistic}}\!\left( t \right)$ respectively is converted based on Table 3 to get $\widehat{NACp}_{min}$ or $NAC{p_{min}}$ respectively. These new values represent the worst possible NACp value with respect to the standard effects and conditions.

7. 7. The last value denoted as $NAC{p_{ref}}$ , respectively $\widehat{NACp}_{ref}$ , is used to determine when the aircraft is no longer jammed. Specifically, it sets the value of the NACp indicator that the aircraft must transmit to no longer be considered impacted by interference. For $NAC{p_{ref}}$ , respectively $\widehat{NACp}_{ref}$ calculation, the most pessimistic estimation of pseudo-range error $\left( {{\sigma _{max}}} \right)$ should be used. The most pessimistic $\left( {{\sigma _{max}}} \right)$ for GPS constellation is ${\sigma _{max}} = 15.6$ meters without SBAS augmentation [34].

   Subsequently, steps 5A. and 5B. are applied in the manner described above and the Equations (6) or (7) is used based on the receiver type. The resulting HFOM value from Equation (6) or (7) is converted into $NAC{p_{ref}}$ or $\widehat{NACp}_{ref}$ respectively using conversion described in Table 3. It should be noted that the recovery of the aircraft’s GNSS receiver can take several minutes after the aircraft leaves the jammed area [11].

2.5 Algorithm for GNSS interference detection

The logic of the algorithm is based on comparing the previously received NACp value at time $t - 1$ with the newly received NACp value at time $t$ . Along with the NACp values, the aircraft’s state at time $t - 1$ , meaning whether the aircraft was jammed or not, plays a key role. As above mentioned, the aircraft is considered not affected by jamming in the first loop of the algorithm when the first ADS-B message is received.

A new loop of the algorithm is executed when a new ADS-B message is received. A comparison of the newly received inputs with the parameters described in Section 2.4 is used to determine the current state of the aircraft. The algorithm workflow in all possible scenarios is described for an aircraft equipped with a GNSS receiver type 3. The detection algorithm runs the same receiver types 1 and 2, only differing in the use of ${\hat \sigma _{max}}$ instead of ${\sigma _{max}}$ , ${\widehat{HFOM}}_{pessimistic}\!\left( t \right)$ instead of $HFO{M_{pessimistic}}\!\left( t \right)$ , and $\widehat{NACp}_{ref}$ instead of $NAC{p_{ref}}$ , as explained in Section 2.4.

2.5.1 Aircraft previously not affected by jamming

In the case, when the aircraft was not affected by jamming in the previous loop, two scenarios can be distinguished. In the first, the updated NACp value increases $NACp\!\left( t \right) \gt NACp\!\left({t - 1} \right)$ . In the second, the updated NACp value drops or remains the same $NACp\!\left( t \right) \le {\rm{\;\;}}NACp\!\left({t - 1} \right)$ .

1. 1. In case $NACp\!\left( t \right) \gt NACp\!\left({t - 1} \right)$ , which means the position error got even smaller than in the last ADS-B message, the methodology automatically marks the aircraft as non-jammed.

2. 2. In case $NACp\!\left( t \right) \le {\rm{\;\;}}NACp\!\left({t - 1} \right)$ , the algorithm estimates whether the aircraft remains non-jammed or becomes jammed. This is conducted by a comparison test between the transmitted $NACp\!\left( t \right)$ and the calculated $NAC{p_{min}}$ , representing the maximal error caused by the current constellation. If the $NACp\!\left( t \right) \gt NAC{p_{min}}$ condition is met, the aircraft remains non-jammed. In the opposite case, when $NACp\!\left( t \right) \le {\rm{\;\;}}NAC{p_{min}}$ , jamming is detected and the aircraft state changes to jammed.

Figure 3 displays the workflow diagram of the detection algorithm. It showcases scenarios for a GNSS receiver type 3 that is not impacted by jamming at the start of the algorithm.

Figure 3. Example of the detection methodology workflow for GNSS receiver type 3 that is not previously jammed.

2.5.2 Aircraft previously affected by jamming

In case the aircraft is evaluated as jammed in the previous step, two scenarios may occur based on whether $NACp\!\left( t \right)\ \geq\ NACp\!\left({t - 1} \right)$ or $NACp\!\left( t \right) \lt NACp\!\left({t - 1} \right)$ .

1. 1. In case $NACp\!\left( t \right)\ \geq\ NACp\!\left({t - 1} \right)$ , it has to be decided whether the aircraft remains in the jammed state. First, the comparison test between the $NACp\!\left( t \right)$ value and the calculated $NAC{p_{min}}$ is carried out. If $NACp\!\left( t \right) \le {\rm{\;\;}}NAC{p_{min}}$ , the aircraft state remains unchanged. If $NACp\!\left( t \right) \gt NAC{p_{min}}$ , another test is conducted to check for the possible end of the interference.

Following $NACp\!\left( t \right) \gt NAC{p_{min}}$ , it is necessary to determine the NACp value at which the aircraft is no longer impacted by interference. This tipping point NACp value is referred to as $NAC{p_{ref}}$ . It represents the most pessimistic value for GPS position error without SBAS augmentation. On the condition that $NACp\ \geq\ NAC{p_{ref}}$ , the end of interference is indicated, and the aircraft status changes to non-jammed. Otherwise, when $NACp \lt NAC{p_{ref}}$ , the aircraft state remains unchanged.

1. 2. In case $NACp\!\left( t \right) \lt NACp\!\left({t - 1} \right)$ , meaning an already jammed aircraft provides an even worse position error, the aircraft automatically remains in the jammed state.

It can be assumed that in the case of interference, the position data transmitted by the aircraft is degraded or completely missing. In order to minimise the effect of this phenomenon, aircraft position can be obtained by using other surveillance technology independent of GNSS. It is also possible to use the aircraft’s last reported position as the area affected by jamming is usually not large enough to cause a significant change in HDOP.

2.5.3 Algorithm outputs

All the outputs of the detection algorithm are shown in Table 4. The first five columns represent the algorithm inputs obtained via ADS-B messages while the rest of the columns represent all computed variables within the algorithm. The $HFO{M_{max}}$ and $HFO{M_{pessimistic}}\!\left( t \right)$ values are considered intermediate calculations, therefore they are not stored as outputs.

Table 4. Example of the methodology calculations and output characterised by Y/N value

In the example in Table 4, the aircraft is equipped with GNSS receiver type 3 and previously not affected by jamming. Based on these factors, the interference detection is performed by the comparison of the $NACp\!\left( t \right)$ and $NAC{p_{min}}$ , $NACp\!\left( t \right)\ \geq\ NAC{p_{min}}$ . The state output of the algorithm represents the non-jammed state of the aircraft as 0 and the jammed state as 1. As can be seen in Table 4, in the second loop the algorithm identified the aircraft as jammed since the condition was not met.

3.0 Validation methodology

The validation of the methodology was performed by leveraging data from a series of experiments with GNSS jamming. The experiments were carried out at LKTB airport on the 12th of November 2020 and the 24th of June 2021, and at LKBE airport on the 24th of June 2022. Experiments were organised by the CTU with the participation of the Czech Telecommunication Office, and the Czech Air navigation service provider which issued NOTAMs informing about the presence of GNSS jamming in the vicinity of the affected airports.

In all experiments, a GNSS jammer model TG5CA was used. Although it is a multiband chirped device only the GNSS L1 channel was active during experiments. The GNSS L1 channel has a bandwidth of 20 MHz centred at 1575.42 MHz and a maximum output power of 36 dBm. The GNSS jamming signal can be described as a rapid sweeping chirp signal with a chirp repetition frequency of 100 kHz. The spectrogram and waterfall diagrams of the jamming signal generated by TG5CA are shown in Fig. 4. To ensure complete ADS-B coverage, several ADS-B receivers were deployed within and in the proximity of the test area.

Figure 4. The spectrogram (left) and waterfall (right) diagrams of the jamming signal generated by the TG5CA jammer. In the spectrogram, the blue represents the authentic GNSS signal, and the green represents the GNSS jamming signal from TG5CA.

Several scenarios were performed during the experiments. These included static jammer placement at the airport, dynamic jammer movements around the airport, and jammer on board the aircraft. The controlled environment of the experiments allowed for close monitoring of the aircraft avionics via video recordings. Together with known jammer location and power, the video recordings were used to determine the exact timeframes when aircraft were impacted by GNSS jamming.

The ADS-B data acquired during experiments were organised into seven sets, each originating from a different aircraft. Including both jammed and non-jammed data, 64,580 NACp values were evaluated in total. All aircraft participating in the experiments were equipped with ADS-B version 2 transponders. The GNSS jamming detection success rate of the methodology is evaluated using the confusion matrix and its metrics as shown in Table 5. The metrics used are defined as follows [37, 38]:

• • True positives (TP) – the number of times the methodology has successfully detected jamming in its presence.

• • True negatives (TN) – true negatives: the number of times the methodology determined that interference was not present, but in fact, jamming was not present.

• • False positives (FP) – the number of times the methodology has determined interference was present, but in fact no interference was present. (Also known as “Type I error.”)

• • False negatives (FN)- the number of times the methodology determined that no interference was present, but in fact interference was present. (Also known as “Type II error.”)

Table 5. Confusion matrix used for validation of the methodology

From the four core metrics that make up the confusion matrix, the following additional metrics can be calculated to describe the success of GNSS interference detection in more detail [38, 39]:

• • True positive rate (TPR) – when interference is actually present, how often a positive result is issued.

\begin{equation*}TPR = \frac{{TP}}{{TP + FN}}\end{equation*}

• • False positive rate (FPR) – when interference is not actually present, how often a positive result is issued.

\begin{equation*}FPR = \frac{{FP}}{{FP + TN}}\end{equation*}

• • Precision (PPV) – when a positive result is issued, how often it is the correct result.

\begin{equation*}PPV = \frac{{TP}}{{TP + FP}}\end{equation*}

• • Accuracy (ACC) – overall, how often the issued result is correct.

\begin{equation*}ACC = \frac{{TP + TN}}{{TP + FN + FP + TN}}\end{equation*}

• • Misclassification rate (MISC) – overall, how often the issued result is false, equivalent to 1-Accuracy, also known as ‘Error Rate’.

\begin{equation*}MISC = \frac{{FP + FN}}{{TP + FN + FP + TN}}\end{equation*}

4.0 Results of methodology validation

Table 6 presents the results from the datasets of six aircraft captured during experiments at LKTB. Among these aircraft, one was not jammed during the experiment, however, there was a fluctuation in the NACp indicator, so it was included in the validation. The other five aircraft were jammed for certain time intervals. The total detection accuracy on the LKTB datasets achieved 98.92%.

Table 6. Results of the methodology validation on data from GNSS RFI experiment at LKTB

In the case of LKBE experiments, only one aircraft, Tecnam P2006T, was involved. Several scenarios were performed when the aircraft was flying in and out of the jamming area at different mutual position configurations between the aircraft and the jammer. The aircraft was equipped with a GNSS receiver with SBAS capability, specifically a Garmin 1000 GNSS receiver. The achieved results are displayed in Table 7. In this case, an accuracy of 82.84% was achieved.

Table 7. Results of the methodology validation on data from GNSS RFI experiment at LKBE

Summing up the results, the total accuracy of the methodology evaluated on different GNSS jamming experiments achieved the value of 98.40%. The misclassification rate remained under 2% over the entire dataset. All the evaluation metrics are showcased in Table 8.

Table 8. Total results of the GNSS RFI detection methodology validation

5.0 Discussion

The proposed methodology demonstrated the feasibility of GPS jamming detection using the ADS-B NACp quality indicator in comparison with the maximum possible change of position accuracy achievable using the current constellation of GPS satellites. The proposed methodology achieved an accuracy of 98.40%, precision of 84.34% and misclassification error of 1.60% in jamming detection. The methodology was verified on data from two deliberate experiments with GNSS jamming performed by the CTU.

A more detailed analysis of the results in Section 4 shows that in some cases the transmitted NACp indicator may be lower than the $NAC{p_{min}}$ value estimated by the methodology despite the aircraft not being affected by jamming. This situation leads to a false positive jamming detection rate of 1.49%. There may be multiple causes that lead to this situation. One of them may be the delay in the recovery of the GNSS receiver positioning accuracy after being jammed. Another possibility that could cause this misdetection is the roll movement of the aircraft [18]. If the aircraft is in a high roll angle turn, some satellites may be obscured by the aircraft’s fuselage. In other words, the calculated HDOP value does not correspond to the satellites actually visible in this situation because some may be obstructed by the wing or other parts of the aircraft.

If NACp indicator values transmitted are higher than the calculated $NAC{p_{min}}$ , the methodology detects no jamming and the aircraft is identified as non-jammed even though the jamming is impacting the aircraft. These cases are indicated as false negatives. The possible cause of such cases may be the aircraft’s high roll, pitch or yaw angles and the relative position of the aircraft and the jammer which may cause part of the aircraft to shield the GPS antenna and ‘protect it’ from the jamming signal.

As aforementioned, several studies have been published on the subject of GNSS jamming detection and jammer location estimation using ADS-B quality indicators. The proposed methodology combines the easily accessible satellite constellation data and the ADS-B position error indicator. This combination avoids false interference detection when the satellite positioning service does not provide the required performance due to the visibility of a small number of satellites or their poor geometry. The presented methodology also enables jamming detection only from short-term drops of the NACp indicator. This makes it possible to determine the presence of an interference signal before the aircraft position information is completely lost. Another benefit in comparison to previous research is that the methodology only requires ADS-B data from a single aircraft to detect interference. The need for only one aircraft being affected to detect jamming allows for much faster interference detection in comparison to previously published detection methodologies [11, 13].

An exact comparison of detection performance metrics like true or false detection rates with other methodologies is not feasible since the data used in the validation differs. Importantly, the presented methodology was validated against a precise benchmark derived from deliberate experiments. Whereas the majority of previously published studies used ‘educated guesses’ of the timeframes when aircraft were affected by jamming derived mostly from pilots’ reports.

The only methodology where an exact comparison can be made is the work by Steiner and Nagy [20] who proposed a discrete mathematical model for GNSS interference detection using ADS-B quality parameters. The comparison is feasible since this model was validated against the same set of ADS-B data. In comparison both methodologies achieved detection accuracy of over 98% and misclassification rate under 2%. The biggest difference can be observed in the precision metric where the discrete mathematical model achieved over 7% better score. On the other hand, the model requires a large amount of training data which is not required for the proposed methodology.

Although spoofing detection is out of the scope of the presented methodology, the evaluation of a spoofed aircraft by the methodology was considered. An asynchronous spoofing attack causing a sudden change of position and time could cause a momentary drop in the quality parameters. This drop could be presented as a short-term jamming event by the methodology since the methodology does not differentiate between jamming and spoofing, as spoofing is out of scope. In terms of spoofing signal altering the navigation message, a change of bits describing the satellite health could lead to arbitrarily low computation of HDOP values, which could cause a decrease in detection sensitivity.

6.0 Conclusion

The ADS-B technology has the potential to detect GNSS RFI by analysing the quality indicators contained in ADS-B messages. GNSS interference presents a serious safety threat to aviation. In the event of a disruption, it is essential to maintain the required navigation performance using supplementary navigational aids. This may result in changes to navigation procedures and increased workload for air traffic controllers and pilots, affecting safety, air traffic capacity and associated costs. Early detection of GNSS RFI can help airspace users prepare for changes in procedures and mitigate the interference impact. Therefore, a new GNSS jamming detection methodology is proposed in this paper that compares the ADS-B NACp quality indicator with HDOP values obtained from the GPS almanac.

The detection methodology involves several steps. First, the HDOP value is calculated based on the GPS satellite almanac data, taking into account all satellites in the aircraft’s line of sight and excluding those that are obscured by the terrain or unhealthy. Next, the HFOM value is derived from the NACp indicator provided by the aircraft. The GNSS jamming detection algorithm estimates the minimum possible $NAC{p_{min}}$ value for the current constellation of satellites and determines whether the transmitted NACp value is lower or higher than the estimated minimum. If the transmitted value is lower, the aircraft is considered jammed. The methodology also explains how to determine the end of jamming detection. The performance of the methodology was evaluated using the data from deliberate GNSS RFI experiments and computing the confusion metric. The methodology achieved a detection accuracy of 98.40%.

The presented methodology provides a useful tool for GNSS jamming detection and is relatively easy to implement into operation. However, the current form has some limitations. In terms of speed, the algorithm is slowed down by downloading data from the API to obtain terrain elevation for a specific latitude and longitude. It is therefore recommended to use a geographical information system tool to identify potential satellite obstructions. It is also possible to avoid calculating aircraft visibility of satellites for cases where the aircraft is flying at higher flight levels. In the future, when dual-frequency multi-constellation GNSS receivers get implemented into aircraft, the methodology can be extended by almanac data from additional constellations.

In summary, the methodology leverages aircraft GNSS receivers as GNSS jamming sensors and can be easily implemented using available ADS-B and GPS almanac data. The interference evaluation algorithm itself is not computationally intensive. Additionally, the methodology does not require any historical interference ADS-B data for model training. Overall, this methodology offers a practical solution for GPS interference detection and achieves satisfactory results. Even better performance could potentially be achieved when combined with other GNSS RFI detection tools or methodologies.