Locating GNSS Interference Sources using ADS-B with Non-linear Least Squares

2.2.1 Free-Space Propagation Loss

The first step in predicting the jammer power that would be received at an aircraft involves modeling the power loss between the transmitter (jammer) and the receiver (aircraft). Assuming a known interference source located at (x, y, z) in an Earth-centered Earth-fixed (ECEF) coordinate frame, and given our assumption of unobstructed signal propagation from the transmitter to the receiver, we can use the Friis equation (Friis, 1946) to predict the jamming power received at a specified aircraft location (x_a, y_a, z_a). This equation is as follows:

$P_r=\frac{G_t{G}_r{\lambda }^2{P}_t}{(4\pi {)}^2{d}^2}$ 1

$\text{where}{d}^2={(x-x_a)}^2+{(y-y_a)}^2+{(z-z_a)}^2$

Here, P_r is the received power in watts, P_t is the transmitted power in watts, λ is the transmitted signal wavelength in meters, d is the distance between the transmitter and receiver in meters, and G_t and G_r are the (unitless) transmit and receive antenna gains, respectively.

For the purposes of this study, expressing power in watts can introduce numerical rounding errors due to the small magnitudes and wide dynamic range of many power measurements, which can range from milliwatts to picowatts. We therefore re-express the Friis equation in dBW by applying 10 log₁₀ to both sides. In addition to mitigating any rounding errors, the logarithmic transformation simplifies the process of taking the derivative of the objective function:

$P_{r_{db}}=G_{t_{db}}+G_{r_{db}}+P_{t_{db}}+20{\log }_{10}\left(\frac{\lambda }{4\pi }\right)-10{\log }_{10}[{(x-x_a)}^2+{(y-y_a)}^2+{(z-z_a)}^2]$ 2

In the objective function in Section 2.3, the parameters x, y, z and P_t become the unknowns that need to be determined.

2.2.2 Antenna Gain Pattern

Because we assume the jamming source is omnidirectional, the transmitted gain (G_t) in Equation (1) is set to one. However, accurately modeling the received gain (G_r) requires a comprehensive representation of the receiver’s gain pattern. To achieve this, our algorithm uses the relative radiation pattern requirements established by the RTCA, as shown in Figure 3 (RTCA SC-159, 2018). These standards outline the minimum operational performance criteria for active airborne GNSS antenna equipment operating in the L1/E1 and L5/E5a frequency bands. We expect that nearly all antennas conform to the RTCA specifications, making them a reliable basis for modeling G_r.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

RTCA requirements for active airborne GNSS antenna equipment

Ideally, GNSS receiver antennas for airborne applications would exhibit zero gain at negative elevation angles. However, achieving this would require an extensive ground plane, which is challenging or impractical in real-world settings. In practice, the gain pattern at negative elevation angles depends on factors such as ground plane size and the resonant patch element design. Our model uses data from Bauregger (2003), who delineated the gain pattern at negative elevation angles for microstrip patch antennas, the most prevalent type of airborne GNSS antenna (Bauregger, 2003). Table 2 presents the relative radiation pattern template taken from Bauregger (2003). This pattern is important because most jamming signals originate from the ground and are therefore most likely to be received by upward antennas at negative elevation angles.

2.2.3 Line-of-Sight Propagation

Due to the design of our objective function for interference localization (Section 2.3), we must also account for scenarios in which the aircraft remains unaffected by the jammer. For example, if the aircraft is not within the line-of-sight of the jammer, the jamming signal will likely be obstructed by the curvature of the Earth’s surface. This consideration is particularly important for powerful jammers, which can impact aircraft at higher altitudes but will not necessarily affect aircraft located at greater horizontal distances from the jammer. Such situations are observed in interference events over regions like the eastern Mediterranean.

Signal propagation in the jammer’s line-of-sight can be modeled as follows:

$\sqrt{{(x_a-x_j)}^2+{(y_a-y_j)}^2+{(z_a-z_j)}^2}\le \sqrt{2\cdot \frac{4}{3}{R}_E{h}_j}+\sqrt{2\cdot \frac{4}{3}{R}_E{h}_a}$ 3

where (x_j, y_j, z_j) is the jammer location in ECEF, (x_a, y_a, z_a) is the aircraft location in ECEF, R_E is the Earth’s radius, h_j is the ground elevation of the jammer, and h_a is the altitude of the aircraft. Figure 4 provides a visual illustration of this line-of-sight signal propagation. If the above condition is met, then the aircraft has a line-of-sight to the jammer, and vice versa. The derivation of Equation (3) is detailed in Appendix A.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Illustration of line-of-sight propagation in ECEF

2.2.4 Relationship between NIC and Received Jamming Power

At this point, our mathematical model depicts the interplay between a jamming source and an aircraft based on signal power. It allows us to predict the received power (${\hat{P}}_r$) at an aircraft given a presumed jammer location and jamming power (yellow-shaded segment in Figure 2). However, the remaining challenge is to estimate the received jamming power (${\tilde{P}}_r$) using information contained in ADS-B messages.

We do this using the NIC parameter, which, as noted in Section 1.1, indicates GNSS signal quality by showing the integrity containment radius within which the aircraft’s current position is guaranteed to reside. This horizontal protection level can be affected by several factors including the standard deviation of measurement noise, satellite geometry, and other considerations like the maximum allowable probabilities for a false alert.

Because RTCA DO-260B (RTCA SC-186, 2011) requires NIC to be ≥ 7 under normal operating conditions, we assume that any instance of NIC < 7 indicates jamming-induced satellite signal tracking loss. We accordingly define the jamming power level corresponding to NIC = 0 as the receiver’s tolerance threshold. If the jamming power exceeds this threshold, the receiver is expected to cease tracking, rendering the reported protection level uncertain.

To determine the appropriate value for tolerable jamming power, we use the derivation by Kaplan and Hegarty (Kaplan & Hegarty, 2017). Figure 5, which is adapted from Kaplan and Hegarty (2017), shows the tolerable jamming power as a function of the C/N₀ tracking threshold under conditions of wideband interference with band-limited white noise. The figure is based on a third-order phase-locked loop model for unaided GNSS receivers. For traditional BPSK-R(1) signals like L1 C/A, a typical C/N₀ tracking threshold is approximately 27 dB-Hz. In this study, we adopt a more conservative threshold of 25 dB-Hz, which is consistent with design standards for aviation-grade receivers and corresponds to a tolerable received jamming power of −115 dBW. We therefore assume that, when NIC = 0, the received jamming power exceeds this −115 dBW threshold.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Tolerable jamming power as a function of the tracking threshold for L1 C/A, L2 CL and L5 Q5 signals under wideband noise interference. Adapted from Kaplan and Hegarty (2017).

Given this assumption for NIC = 0, the next step is to analyze the empirical relationship between received jamming power and NIC for NIC ∈ [1,6]. To explore this relationship, we developed a testbed that emulates the RFI environment experienced by airborne GNSS receivers (Chen et al., 2023). In our setup, live-sky GPS signals are captured and re-radiated within an anechoic chamber using an aviation-grade antenna. These re-radiated signals are received by multiple GNSS receivers positioned outside the chamber. Each receiver collects raw measurements, which are used to compute the horizontal protection level at each time step using a standard Receiver Autonomous Integrity Monitoring (RAIM) algorithm (Blanch et al., 2015). NIC values are then derived from the computed horizontal protection level according to Table 1.

Jamming signals are introduced into the system using a programmable transmitter connected to a separate antenna within the chamber. The transmitted jamming power is controlled via the transmitter gain, while the actual received power P_r is measured at the input of a power meter. We tested two jamming types: wideband noise and dual-tone continuous wave. For each jamming scenario, we collected approximately 30,000 observations using each of three commercial GNSS receivers: the Septentrio Mosaic-X5, Trimble BX940, and u-blox ZED-F9P.

Figure 6 presents scatter plots showing the relationship between NIC values and received jamming power across the different receiver brands and jamming signal types. A small amount of random jitter was added to the plotted points to minimize overlap. Each circle represents a unique NIC value observed at a given P_r level. Sparse clusters, reflected by incomplete or faint circles, correspond to infrequently observed value combinations.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

NIC values given different received jamming power from continuous wave (CW; left) and wideband noise (right) interference at three different GNSS receivers (rows)

As shown in the right-hand side of Figure 6, all three receivers exhibit broadly similar susceptibility and tolerance to wideband noise jamming. In contrast, their responses to continuous wave jamming differ significantly: the Trimble receiver shows little resilience, while the u-blox receiver demonstrates substantial resistance, especially at lower power levels (below −130 dBW). These differences highlight that, while general trends exist, receiver-specific responses introduce notable variability. Because of this variability among receivers, it is difficult to deterministically associate specific NIC values with specific received power levels. Nevertheless, for NIC values between 1 and 6, the data indicate that the corresponding received jamming power remains within a narrow range of 3 to 5 dBW across all receivers.

Figure 7 shows how NIC values also depend on the number of satellites tracked by the receiver. Our testbed results indicate that NIC ∈ [1,6] typically occurs when exactly five satellites are used, while NIC = 0 is reported when only four or fewer satellites are available. However, even when the satellite count remains constant (e.g., at five), environmental factors such as aircraft orientation, body shielding, and interference directionality can still affect NIC values. This variability highlights how NIC is influenced by both satellite geometry and signal quality.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

NIC values versus the number of satellites tracked by different receivers (rows) experiencing continuous wave (left) and wideband noise (right) interference

Given these various sources of complexity for NIC ∈ [1,6], we propose using NIC as a binary indicator of jamming severity. Specifically, we adopt a C/N₀ tracking threshold of 25 dB-Hz for unaided receivers, consistent with industry norms (Kaplan & Hegarty, 2017, Table 9.7), which corresponds to a tolerable received jamming power of −115 dBW for L1 C/A signals under wideband interference. We accordingly define:

• • $\text{NIC}=0\to {\tilde{P}}_r\ge -115dBW$ (severe jamming; tracking lost)

• • $\text{NIC}\in \left[\text{1,6}\right]\to {\tilde{P}}_r<-115dBW$, typically within a 3−5 dBW range (moderate jamming; tracking degraded)

While this binary mapping is empirically derived and less precise than direct C/N₀-based measurements, the trend is sufficiently robust to infer interference. Any unmodeled offsets, such as variation in receiver design, antenna gain, or signal obstruction, can be absorbed into the estimated transmit power P_t. In our framework, P_t represents an effective or relative jamming source strength and can therefore be interpreted as the true transmit power plus an unknown receiver-dependent bias.

2.3.1 Non-linear Least Squares Problem

The first step in our solution involves formulating an objective function. Given a subset of ADS-B data, assume there are m data points for a single flight. Each data point includes the aircraft’s location in ECEF coordinates and the estimated received jamming power (${\tilde{P}}_r$) inferred from the NIC value: $(x_a,y_a,z_a,{\tilde{P}}_{ra})\in {ℝ}^m$. The unknown parameters that must be estimated are the jammer’s location and transmitted jamming power: (x, y, z, P_t). Starting with an initial guess for the jammer’s location and power, we can use Equation (2) to predict the received jamming power at each aircraft location. Candidate jammer locations are typically distributed across a coarse grid covering a large bounding box that contains all observed low-NIC points. For each grid point, we also evaluate possible transmit powers P_t ranging from 0.01 watts to 1000 watts, which should cover interference sources ranging from small portable GNSS jammers to high-power military-grade systems. This broad range of candidate jammer powers ensures a robust initial search in our optimization procedure.

The objective function f is defined as the weighted sum of squared residuals between the estimated (${\tilde{P}}_r$) and predicted (${\hat{P}}_r$) received jamming power, where the residuals are derived from the differences specified in Table 3:

$\mathrm{minimize}f\left(x\right)={({\stackrel{\sim }{P}}_r-{\hat{P}}_r)}^{T}W({\stackrel{\sim }{P}}_r-{\hat{P}}_r)$ 5

where W is a weight matrix that allows us to design the mechanisms behind the model:

$W=R^{-1}={\left[\begin{array}{ccccc}{\sigma }_2^2& {\rho }_{12}{\sigma }_1{\sigma }_2& & {\rho }_{1(m-1)}{\sigma }_1{\sigma }_{(m-1)}& {\rho }_{1m}{\sigma }_1{\sigma }_m\\ {\rho }_{21}{\sigma }_2{\sigma }_1& {\sigma }_2^2& & {\rho }_{2(m-1)}{\sigma }_2{\sigma }_{(m-1)}& {\rho }_{2m}{\sigma }_2{\sigma }_m\\ {\rho }_{31}{\sigma }_3{\sigma }_1& {\rho }_{32}{\sigma }_3{\sigma }_2& & {\rho }_{3(m-1)}{\sigma }_3{\sigma }_{(m-1)}& {\rho }_{3m}{\sigma }_3{\sigma }_m\\ ⋮& ⋮& \cdots & ⋮& ⋮\\ {\rho }_{(m-1)1}{\sigma }_{(m-1)}{\sigma }_1& {\rho }_{(m-1)2}{\sigma }_{(m-1)}{\sigma }_2& & {\sigma }_{m-1}^2& {\rho }_{(m-1)m}{\sigma }_{(m-1)}{\sigma }_m\\ {\rho }_{m1}{\sigma }_m{\sigma }_1& {\rho }_{m2}{\sigma }_m{\sigma }_2& & {\rho }_{m(m-1)}{\sigma }_m{\sigma }_{m-1}& {\sigma }_m^2\end{array}\right]}^{-1}$

In our algorithm, the weight matrix is the inverse of the covariance matrix (R). The diagonal elements (${\sigma }_i^2$), which represent the variance of the i-th measurement, encapsulate two key factors. First, they capture the inherent uncertainty in inferring jamming power (${\tilde{P}}_r$) from NIC, as outlined in Section 2.2. Second, they reflect the informativeness of each measurement:

1. Points affected by interference are assigned more weight than unaffected points. This decision is driven by the challenge of data imbalance: when investigating interference, the selected airspace for analysis typically exceeds the true interference region to ensure thorough coverage, resulting in more unaffected points than affected ones.

2. Points corresponding to the initial drop in NIC are also given greater emphasis. These points help identify the boundaries of the interference region, as the initial drop in NIC occurs as aircraft approach and enter the affected area. In contrast, as aircraft leave the affected area, the receiver may take time to recover from interference, causing a more delayed return to accurate position solutions and normal NIC values. Because this delayed response can lead to low NIC values in regions unaffected by interference, we assign less weight to points leaving the affected region.

The off-diagonal elements (ρ_ijσ_iσ_j) capture autocorrelation between measurements from the same aircraft, as data points that are close to each other in time are often highly similar. Here, we define a change point as when the interference event begins or ends, or the times when the aircraft enters or leaves the affected region. The correlation between two temporally consecutive points i and j that both occur either before or after this change point is modeled as ${\rho }_{ij}=e\frac{-|t_i-t_j|}{\tau }$, where τ is a hyper-parameter determining the duration over which two points remain highly correlated. If the time difference between two points exceeds τ, we treat them as independent. We chose a τ value of 20 sec because interference events in KDEN occur within a relatively small area. For more powerful or longer-lasting interference events, τ may increase substantially, as more measurements will be redundant over longer durations.

2.3.2 Gauss Newton

Given the non-linear objective function described above, our next step is to determine the global minimum point using the Gauss-Newton method. The pseudocode is as follows:

In Line 3 of the above pseudocode, we linearize the objective function (Equation (5)) near the local state at each iteration using a first-order approximation:

$f\left(x\right)\approx f\left(x^{\left(k\right)}\right)+Df\left(x^{\left(k\right)}\right)(x-x^{\left(k\right)})=A^{\left(k\right)}x-b^{\left(k\right)}$ 6

where Df(x^(k)) = the Jacobian matrix, with entries ${\left(Df\right)}_{ij}=\frac{d{f}_i}{d{x}_j}$

$\begin{array}{cc}& k=\text{the}{k}_{\text{th}}\text{iteration}\\ & A^{\left(k\right)}=Df\left(x^{\left(k\right)}\right)\\ & b^{\left(k\right)}=Df\left(x^{\left(k\right)}\right)x^{\left(k\right)}-f\left(x^{\left(k\right)}\right)\end{array}$

The matrices A^(k) and b^(k) are introduced to simplify the first-order approximation equation and represent the solution of the linearized least squares problem in a more familiar form. From Equation (2), each element of the Jacobian matrix Df(x^(k)) can be calculated by finding the gradient of received power with respect to each variable x, y, z, and P_t:

$Df\left(x^{\left(k\right)}\right)=10\left[\frac{2(x_j^{\left(k\right)}-X_a)}{\ln \left(10\right)d^2}\frac{2(y_j^{\left(k\right)}-Y_a)}{\ln \left(10\right)d^2}\frac{2(z_j^{\left(k\right)}-Z_a)}{\ln \left(10\right)d^2}\frac{-1}{\ln \left(10\right)P_t^{\left(k\right)}}\right]$ 7

where $f\left(x^{\left(k\right)}\right)=P_r-(G_{t_{db}}+G_{r_{db}}+10{\log }_{10}{P}_t^{\left(k\right)}+20{\log }_{10}\left(\frac{\lambda }{4\pi }\right)-10{\log }_{10}{d}^2)$

$d^2={(x_j^{\left(k\right)}-x_a)}^2+{(y_j^{\left(k\right)}-y_a)}^2+{(z_j^{\left(k\right)}-z_a)}^2$

In Line 4 of the above pseudocode, we solve the linearized objective function A^(k)x − b^(k). Because we have more measurements than unknown parameters, the matrix A^(k) is skinny and full rank. The least-squares solution can be obtained by setting the gradient with respect to x to zero, yielding:

$X^{(k+1)}={\left(A^{{\left(k\right)}^T}W{A}^{\left(k\right)}\right)}^{-1}{A}^{{\left(k\right)}^T}W{b}^{\left(k\right)}$ 8

2.3.3 Error Bound

Given the assumptions outlined in Section 2.1, our algorithm may encounter a degradation zone in which it fails to accurately locate the jamming source. To identify the boundaries of this region, we calculate a 95% confidence interval by computing the covariance matrix of the estimated jammer state. This matrix is computed using Equation (9), as described in Appendix B:

$Var\left(\hat{X}\right)=E\left({\hat{X}}^2\right)-E{\left(\hat{X}\right)}^2={\left(A^{T}WA\right)}^{-1}$ 9

Using Equation (9), the confidence interval for the estimated jammer state becomes:

${\hat{x}}_i\pm c\sqrt{\mathrm{var}\left({\hat{x}}_i\right)}$ 10

where $\mathrm{var}\left({\hat{x}}_i\right)$ denotes the i_th diagonal element of (A^TWA)⁻¹ in Equation (9), and c represents the confidence factor associated with the chosen confidence level. The derivation of Equation (10) is detailed in Appendix C. In this study, we estimate the 95% confidence interval, for which the corresponding c-value is 1.96 given a standard normal distribution.