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 [Reference Ma and Zhou28]

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 [Reference Ma and Zhou28]

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 [Reference Powell30]. 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 [Reference Harriman31], 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 [Reference Elevation32] 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 [Reference Elevation32] 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 [Reference Jonas and Vitan11].

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.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.