Optimal INS Monitor for GNSS Spoofer Tracking Error Detection

In this article

target receiver's tracking loops to transition from the authentic to the spoofed signals and to then slowly drag the user along a false trajectory.A spoofer can also opt to jam the GNSS signals before transmitting the spoofed signal.A more sophisticated spoofer could send two spoofed signals, where one of the spoofed signals would be in opposite phase to the authentic signals, causing cancellation.In this process, known as nulling, the target receiver would not even be able to detect multiple correlation peaks.Ideally, spoofing detection techniques should protect against all possible scenarios.
Our group introduced the first stochastic description and quantification of the performance of an IMU-based GNSS spoofing monitor against worst-case "faults" (i.e., spoofing inputs over time) (Khanafseh et al., 2014;Tanil, Khanafseh, et al., 2018;Tanil et al., 2015aTanil et al., , 2015b;;Tanil, Khanafseh, et al., 2016;Tanil et al., 2017).We specifically investigated anti-spoofing solutions utilizing IMUs because essentially all modern aircraft are equipped with IMUs, thereby requiring minimal additional cost or system modification.An IMU is naturally immune to external interference, which makes it an excellent resource for ensuring navigation continuity.Additionally, when used in the navigation solution in various integration schemes with GNSS (uncoupled or loosely, tightly, or ultra-tightly coupled), the INS provides the redundancy needed to resist spoofing attacks.In our prior work (Tanil, Khanafseh, et al., 2018;Tanil et al., 2017), we developed a chi-squared innovation sequence detector to monitor the accumulated time history of normalized Kalman filter (KF) innovations.The two main advantages of this cumulative innovation (CI) sequence monitor are that innovations are already available in the KF, such that little additional computation is required for the monitor implementation, and that it provides detection capability against slowly growing faults.We evaluated the performance of the CI monitor against worst-case GNSS fault profiles both analytically and experimentally (Tanil, Jimenez, et al., 2018;Tanil, Khanafseh, et al., 2018).The worst-case fault here represents a spoofed GNSS signal profile that maximizes integrity risk.We also analyzed the sensitivity of the CI monitor against error modeling uncertainties in the INS/GNSS KF structure (Kujur et al., 2019).
However, post-detection recovery has not been addressed in previous workthe difficulty being that the KF is already corrupted once spoofing is detected.Moreover, previous performance evaluations assumed that the CI monitor started at spoofing onset and operated without a defined run time.Finally, the CI monitor did not provide the means to produce a protection level-i.e., a position-domain containment boundary corresponding to the maximum acceptable level of integrity risk.To address these critical limitations, in this paper, we introduce a new type of CI monitor: a cumulative position-domain innovation (CPI) monitor that detects spoofing by accumulating the target position tracking error embedded in the spoofer's signal.We also present a complementary solution separation (SS) concept to produce protection levels and provide a means for post-detection recovery.
Section 2 provides useful background information on the original CI monitor.Section 3 introduces the new CPI monitor.An example application with quantitative results is presented in Section 4. In Section 5, we present the complementary SS monitor.Finally, we summarize our work in Section 6.We also provide relevant derivations in the appendices.

KF STATE MODEL
We consider a vehicle employing INS and GNSS sensors integrated with a KF to estimate its position, velocity, and attitude.The dynamics of the INS/GNSS system, augmented as needed with sensor error state dynamics, are linearized to obtain the process model utilized in the KF: where x k is the state vector, Φ Φ k is the state transition matrix, Γ Γ w k is the process noise model matrix, and w k is the additive white process noise with a respective covariance matrix Q k .The measurement model is as follows: where H k is the observation matrix and ν ν k is the measurement noise with a respective covariance matrix V k .The innovation vector γ γ k at time epoch k is defined as follows: where x is the state vector estimate prior to the measurement update at time epoch k.The CI detector is a chi-squared monitor that utilizes the cumulative normalized innovations from a KF as the test statistic and compares this test statistic against a threshold.A cumulative test statistic q N at time epoch N is the sum of squares of the normalized innovation vectors over time, given as follows: where is the innovation vector covariance matrix at time epoch k and P k is the estimate error covariance matrix prior to the measurement update at time epoch k.The monitor simply checks whether the test statistic q N is smaller than a predefined threshold T N 2 .For a given false alarm requirement under fault-free conditions, the threshold T N is determined from the inverse chi-squared cumulative distribution function (CDF) with N degrees of freedom.The monitor produces an alarm if q T N N > .One of the major limitations of the CI monitor is that it does not offer a means for post-detection recovery because there is no fault-free source upon which to rely.Another unaddressed issue concerns the monitor start and run times.In all prior performance evaluations, it was assumed that the monitor start time (conveniently) coincided with the spoofing onset time, and no guidance was provided to determine how long the monitor should run before resetting.In the next section, we introduce the CPI monitor and show that it directly addresses the run time issue and also provides superior detection performance relative to the CI monitor for any run time.The start time and recovery are addressed afterward.

Spoofer's Tracking Error
The initial objective of a smart spoofer is to cause the target receiver to lose lock of the authentic GNSS signals and lock onto the counterfeit signals without being detected.This initial lock transition is crucial because any abrupt changes during the process would be easily detectable, sabotaging the spoofer's plan.The best way for a spoofer to cause the target to switch the counterfeit signal without being detected is to initially replicate the authentic signals at lower power and then slowly increase the power to cause the target receiver to transition.If successful, the spoofer would then attempt to inject small but accumulating position and/or time offsets in an attempt to slowly pull the target away along the desired spoofed trajectory while remaining undetected.
To deliver a replica of the authentic signal to the target, the spoofer would need to know the position of the target's GNSS antenna during the attempted takeover.Consequently, the spoofer would need to track the target in real time.Any tracking errors would ultimately be embedded in the replica signal and appear to the target as additional "noise" in the received GNSS signals.As tracking errors are inevitable, detecting the presence of such unusual noise would expose active spoofing.This is true even during the initial takeover phase when the spoofer has not yet injected any additional offsets to the spoofed signals.In the remainder of this paper, we address this case specifically, as any additional offsets to the spoofing profile are irrelevant to the monitor.The spoofer's uncertainty of the target position will simply be referred to as "tracking error" in the remainder of the paper.In Appendix C, we show that the tracking error appears in the spoofed measurements as an additive term, and the measurement equation can be represented as follows: where z k s is the spoofed measurement vector, z k is the "true" (unspoofed) measurement vector, H k is the observation matrix that maps the tracking error to the range domain, and ν ν k t is an n ×1 column vector of tracking error, with n being the number of states.The superscript s indicates "spoofed." We can observe these tracking errors in the KF innovation vector because it extracts the difference between the GNSS measurements and the predicted measurements of the process model from the INS.We will see later that even small deviations over time are readily observable when the INS is used because of the precision of GNSS carrier phase measurements.
Because the innovation vector γ γ k in Equation ( 3) is constructed from the entire measurement vector z k and the entire state vector x k , the effect of the position tracking errors, our proxy for spoofing, will not be the only contributor to γ γ k .A more direct way to observe the effect of the tracking errors would be to project the innovation vector into the position domain.

Neyman-Pearson Optimal and Sufficient Test Statistic
We desire to find a test statistic to maximize the probability of detection for a given false alarm rate.The Neyman-Pearson lemma (Neyman & Pearson, 1933) provides the solution to this problem: given two mutually exclusive hypotheses H 0 and H 1 , which for some observation x have conditional probability densities p x H ( | ), the likelihood ratio is the optimal test statistic.The likelihood ratio is defined as follows: Our two hypotheses are the fault-free case, which has no tracking error, and the spoofed case, in which the tracking error appears in the innovation vector.As a starting point, we model the tracking error ν k t as white Gaussian noise (WGN) distributed as  (0, ) 2 σ t , where σ t 2 is the unknown variance of the tracking error.We consider the innovation vector distributions with and without the tracking error for our likelihood ratio test.The fault-free innovation vector distribution at any time k is given as H k k 0 (0, ) : γ γ   S , whereas the spoofed innovation vector distribution, as shown in Appendix C, is given as , where R uu T t V and u is a unit vector in an arbitrary spatial direction in which the tracking error exists.The innovation vectors are mutually independent over time under H 0 .In principle, the vectors can also be modeled as independent under H 1 given that a "high-quality" IMU is being utilized because x k in Equation ( 2) would be relatively unaffected by the spoofer's tracking error if a high-quality IMU were used.(Note: We refer to IMUs with specifications equivalent to or better than the navigation-grade IMU listed in Appendix F as "high-quality" and IMUs with specifications equivalent to or lesser than the automotive-grade IMU listed in Appendix F as "low-quality/grade.")We accept this assumption as true in the following derivation of the optimal test statistic, and in Appendix D, we show that this assumption is valid even for small tracking errors and "low-grade" IMUs.
Using Equation ( 6) and given an arbitrary threshold λ( ) N for independent innovations J J k k N ( 1, ..., ) , the likelihood ratio test can be expressed as follows (Kay, 1998): where: captures the terms that do not depend on J J k k N ( 1, ..., ) .Taking the natural log of both sides and collecting all of the constant terms on the right, we obtain the following: Moreover, because the matrices S k and S H RH + are non-singular, we can write the following (Miller, 1981): Substituting Equation (10) into Equation ( 9), we obtain the following relation: where the constant on the right side is as follows: For short periods of time, matrices H and S can be assumed to be time-invariant; thus, the term tr H RH S 1 can be moved to the right side of Equation ( 11).Equation ( 11) then reduces as follows: where: Recalling that R uu T t V , we substitute this term into Equation ( 13) to obtain the following: Rearranging the terms of Equation ( 15) and noting that S −1 is symmetric, we can now write the following relation: where we have defined C N C N t 4 3 ( ) ( )/ V .Finally, we introduce the following scalar projection of the innovation vector: which is a weighted projection of the innovation vector into the position-domain direction u, the tracking error direction under consideration.Therefore, the optimal test statistic is given as follows: which we call the CPI.We note that q N u does not require a knowledge of σ t .

Position-Domain Innovation
Under spoof-free conditions, the scalar position-domain innovation in Equation ( 17) is distributed as follows: To simplify the notation, we define the variance as follows: When the spoofer sends a signal mimicking the authentic signal with some inherent additive tracking error, the latter can be observed in the position-domain innovation.In Appendix C, we show that the innovation for the spoofed case with tracking error for an arbitrary spatial direction u is given by the following: Substituting Equation ( 21) for γ γ k in Equation ( 17), we obtain the spoofed position-domain innovation: The mean of γ k us is again zero, but the variance is now as follows: This can be written more compactly using Equation ( 20): Thus, under spoofed conditions, the position-domain innovation has the following distribution: For notational simplicity, we also define the following: The tracking error affects the position-domain innovation by increasing its variance without changing the mean (which remains zero).Now, we can write the position-domain innovation before and after spoofing, respectively, as follows: When normalized by V J k u , the position-domain innovation is distributed in the unspoofed case as follows: In the spoofed case, the position-domain innovation is distributed as follows:

Sum of Squares of Normalized Position-Domain Innovations
The square of a scalar normal random variable distributed as  (0, ) 2 σ follows the Gamma distribution 1 2 2 , 2V (Mathai & Provost, 1992).The Gamma distribu- tion ( , ) D T is defined by its shape parameter α and scale parameter θ .Therefore, we can write the distribution of the square of the normalized position-domain innovation in the unspoofed case as follows: In the spoofed case, we obtain the following: For a period of accumulation N, following the results in Sections 3.2 and 3.3, we define our optimal CPI test statistic (in the unspoofed case) as follows: We note that the sum of independent Gamma-distributed random variables is also Gamma-distributed (Mathai & Provost, 1992).We assume that the change in variance of the position-domain innovations is negligible over the accumulation period, i.e., V V test statistic in the unspoofed case q N is then Gamma-distributed as follows: For a given probability of false alarm P FA , we can determine the threshold as follows: where F * 1 is the inverse CDF of the Gamma distribution.In the spoofed case, with the tracking error embedded in the test statistic, we have the following: By defining the ratio : , we can rewrite the above equation as follows: From Equations ( 34) and ( 37), we can see that tracking error causes the scale parameter of the test statistic distribution to change but does not affect the shape parameter, which remains the same as in the unspoofed distribution.The probability of missed detection (i.e., of not detecting the tracking error) is as follows: where T N is the threshold, as defined in Equation ( 35), γ( , ) a b is the lower incomplete Gamma function: and Γ( ) a is the Gamma function: * J ( ) ( , ) a a f .

Gaussian Approximation
For a Gamma distribution with shape parameter α and scale parameter θ , the mean and variance are DT and DT 2 , respectively.To determine how the tracking error helps with detection, we take the large-N approximation for a Gamma distribution: Thus, the distributions of q N u and q N us can be approximated for large N, respectively, as follows: We now introduce a modified test statistic: The unspoofed distribution of this new test statistic is as follows: The spoofed test statistic distribution is given as the following: It is clear from the preceding two equations that the accumulated tracking error causes the mean of the test statistic distribution to grow as N increases while the variance of the distribution remains independent of N.

Tracking Error as Colored Noise
In the development thus far, for clarity, we have assumed that the tracking error is WGN.However, it is possible that the tracking sensor error could be time-correlated or even that the spoofer would choose to filter the tracking sensor output to try to smooth out the errors to lower the possibility of detection.
To analyze the performance of the CPI monitor, we need to determine the distribution of the test statistic in the case of colored tracking error.We model the tracking error as a zero-mean first-order Gauss-Markov random process (GMRP) with time constant τ t and variance σ t 2 .To the spoofer's advantage, we ignore any time delays that might be incurred in the filtering process.
Although the CPI monitor operates sequentially in time, to analyze the monitor performance, it is more convenient to use a batch method for quadratic forms of random variables, as described by Mathai & Provost (1992).Let X be an N ×1 random vector distributed as X   (0, ) Σ and A be an arbitrary N N × symmetric matrix.Then, the quadratic form Q X X AX T ( ) = can be expressed as a linear combination of independent central chi-squared variables Y: where λ λ λ λ 1 2 3 , , , ..., N are the eigenvalues of Σ Σ , and P is the orthonormal matrix of eigenvectors of Σ Σ 1/2 1/2 A .We can now express the cumulative position-domain innovation from time 1 to N in terms of X, with Σ defined as follows: and matrix A defined as follows: making Q X ( ) the sum of the squares of the normalized position-domain innovations, which is our CPI test statistic defined in Equation (36).For authentic signals, only the first term in Equation (47) exists because the innovations under normal conditions are white (and V J ' us 2 0 ).In the case of spoofing, Q X ( ) includes colored noise, and in light of Equation ( 46), it will have a generalized chi-squared distribution, whose CDF can be evaluated using the method described by Das & Geisler (2016).

QUANTITATIVE RESULTS
We first performed en route simulations of an aircraft utilizing a navigation-grade IMU (Appendix F) and single-frequency Global Positioning System (GPS) measurements.The aircraft level flight was simulated to start from 41°50'10″ N, 87°37'30″ W with a cruising speed of 454 knots at an altitude of 40,000 ft.The GNSS measurements were generated using the GPS constellation (SC-159, 2020).The error models for the GPS measurements are defined in Appendix A, and the IMUs are listed in Appendix F. The spoofer tracks the aircraft, with errors, to generate and broadcast counterfeit spoofed signals.The simulated spoofed signals are exact replicas of the authentic signals with either additive zero-mean WGN or zero-mean colored noise modeled as a first-order GMRP.In this example, we used a maximum monitor run time of 180 s, a GPS update frequency of 2 Hz (i.e., N ≤ 360), and a false alarm requirement of 10 5 − .Although tracking errors would exist in all three spatial dimensions, we conservatively assumed that the errors are only in the vertical direction and ignore any errors in the other two directions in the analysis.Figure 1 (left) shows the results of the simulation run with the tracking error modeled as WGN with standard deviation V t 2 cm.The superior performance of the CPI monitor relative to the CI monitor is plainly evident in this simple example.Figure 1 (right) also illustrates the performance of the CPI monitor for WGN tracking errors with larger magnitudes.
The results of the direct simulations are instructive; however, such simulations are an unfeasible means for computing missed detection probabilities.Instead, we turn to Equation (38) to quantify the analytical relationship between monitor run time N, WGN tracking error standard deviation σ t , and probability of missed detection P MD .The results, for the same en route aircraft scenario, in Figure 2 show that the missed detection probability decreases with increasing run time as more tracking errors are accumulated over time and shift the mean of the test statistic distribution, as predicted in Equation ( 42).Similarly, an increasing tracking error magnitude (standard deviation) contributes to a shift in the mean of the test statistic, which reduces P MD .Figure 3 shows two-dimensional cuts of Figure 2 for different tracking error variances and monitor run times.For comparison, we computed the missed detection probabilities for different grades of IMUs. Figure 4 shows the variation in the probability of missed detection as a function of monitor run time N for different IMU grades.The WGN tracking error standard deviation for this example was 10 cm.The figure shows that, even with a relatively low-grade (automotive) IMU, the degradation in monitor performance is negligible.
The sensitivity of the CPI monitor to tracking error can be attributed to the precision of GPS carrier phase measurements and the IMU's accelerometer quality, specifically the velocity random walk (VRW).If tracking errors are smaller in magnitude than the relative position errors obtainable from time-differenced carrier phase measurements, the monitor threshold will be too loose to reliably detect a tracking error.Similarly, if the VRW is large, then the position estimate drift prior to the GPS measurement update will be greater than the impact of the tracking error on the measurement itself.Figure 5 illustrates this case by comparing the performance of the CPI monitor in a low-tracking-error scenario, V t 2 cm, for different IMU grades against a "lousy" IMU with a VRW 100 times greater than that of an automotive-grade IMU.The reason that the automotive-grade IMU performs so well is that the accelerometer noise parameters are comparable to those of higher-grade IMUs.The latter are designed primarily to improve gyroscope bias stability, which is a critical performance parameter for controlling drift over longer time periods-significantly longer than those relevant to the spoofing monitor.
We now turn our attention to quantifying the performance of the monitor in the presence of time-correlated (colored) tracking error, coming either from the sensor itself or as a result of filtering of white sensor output by the spoofer.Figure 7 presents example tracking errors over time with the same standard deviation, but different time constants, including WGN.Using the method described in Section 3.6, we obtain the performance results in Figure 6, which show that time-correlated tracking error would help the spoofer remain undetected for longer, with larger filter time constants leading to a higher probability of missed detection.
We also evaluated the monitor performance for different start times during a standard 24-h period (SC-159, 2020) to evaluate the monitor sensitivity to different FIGURE 4 Probability of missed detection P MD versus monitor run time (N) for different IMU grades, with V t 10 cm FIGURE 5 CPI monitor performance with V t 2 cm for navigation-, tactical-, and automotive-grade IMUs (left) and automotive versus "lousy" IMUs (right) satellite geometries.Figure 8 shows the probability of missed detection versus monitor run time for tracking error modeled as WGN with a standard deviation of 10 cm, with 24 curves representing the performance results at 1-h intervals.It is evident that the monitor is effective over the entire day, although there is obviously some variation in performance with satellite geometry.
These results demonstrate that if a CPI monitor is implemented, the conjecture that INS-based spoofing detection is vulnerable to slowly deviating counterfeit signals can be largely dismissed.The results show that the spoofer's target tracking error should be easily detectable as long as the duration of spoofing lasts exceeds a minimum time defined by the variance and time constant of the tracking error.For the example application considered, Figure 6 shows that even for tracking error standard deviations as small as 10 cm and correlation time constants as large as 40 s, the probability of missed detection is negligible after less than 1 min.We have also experimentally validated the INS monitor, as reported by Kujur et al. (2023), using real GNSS and INS data.

FIGURE 7 Illustration of tracking error versus time for white and colored noise processes
The CPI monitor has the notable limitation that P MD -or equivalently, a position-domain protection level-is difficult, if not impossible, to establish in real time.Thus, the CPI monitor performance must be evaluated in advance by simulation over a wide variety of satellite geometries and tracking error characteristicsfor example, by following a more extensive version of the evaluation presented in this paper for the specific GNSS/INS integration and application considered.In the following section, we introduce a complementary SS monitor concept to ensure the detection of faster-onset spoofing, provide the means to produce protection levels in real time, and allow the navigation function to default to INS-only coasting to preserve continuity (at lower accuracy levels as the INS drifts over time).

SS MONITOR WITH SEQUENTIAL WINDOWS
The SS monitor measures the difference in position solutions between the potentially spoofed integrated INS/GNSS KF output ˆk KF x and an ostensibly clean "coasting" INS-only solution x C k .The test statistic at any time k is defined as follows: The variance for the test statistic is given by the following:  continually increase to maintain a constant false alert probability.Thus, as the coasting window length increases, the SS spoofing detection performance will degrade, allowing slowly growing faults to go undetected.While this is obviously not a desirable performance characteristic, it is precisely over these increasing time windows that the CPI monitor will perform the best.Therefore, given the complementary nature of the two monitors, it is natural to consider running them in parallel, where the CPI monitor would detect slowly growing spoofing profiles and the SS would detect more rapidly growing faults.
For the CPI monitor, a minimum window length N min is required to achieve a desired P MD + .Here, the subscript "+" designates the missed detection (integrity risk) requirement allocation for spoofing attacks lasting longer than N min .The CPI monitor threshold would be obtained by using the false alert (continuity risk) allocation for the monitor, P FA + .For shorter attack windows, the SS monitor would provide the means for detection with the allocated missed detection and false alert probabilities, P MD − and P FA − .Obviously, the sum of the two missed detection probabilities must be smaller than the total integrity risk allocated to undetected spoofing.The same is true for the sum of the two false alert probabilities relative to the total continuity risk allocated to spoofing monitoring.
The SS monitor can provide a protection level for spoofing in the spatial direction u as follows: where: and the multipliers k FA − and k MD − are determined from the SS false alert and missed detection requirement probability allocations.The SS detection threshold is k FA SS k V .The protection level will increase with the SS window length because of the growth in σ C k over time.The minimum run length N min for the CPI monitor sets the upper limit on the SS monitor time window.The combined monitor system is implemented using consecutive fixed-length windows of length N min (Figure 9) with an SS monitor initialized at the beginning of each window and terminated at the end of each window.A new SS monitor window is opened at each new GNSS measurement epoch and closed N min epochs later.Figure 9 illustrates this idea, showing the (conceptual) protection level.Because each window has its own INS-only solution (from the SS monitor), if spoofing is detected in a particular window at time t d , the final INS/GNSS KF solution and associated PL from a prior window closed (without detection) any time before t d can be safely used to initialize subsequent fault-free coasting.
At any given time point, a set number of monitor windows will be running, and the false alert requirement allocation for each monitor can be equally divided among these monitor windows to determine the thresholds for each.This approach is conservative, as the test statistics for the different monitors will be correlated, but it is easy to implement.
The maximum protection level produced by any SS monitor will occur at the end of the window, k N min = , because the INS coasting errors, σ C k , will be the largest at that point.
The value of N min is determined by the CPI monitor performance.For the same example GPS/INS implementation and satellite geometry considered in Figure 6, Figure 10 shows the values of N min for different tracking error magnitudes with the missed detection requirement set to 10 7 − .As in Figure 6, it is clear that colored tracking error takes longer to detect than white noise.Using these values of N min , we can now determine the maximum protection levels.Figure 11 shows the results for a false alarm requirement of 10 5 − and missed detection requirement of 10 6 − , based on Equation ( 51).Again, it is clear that the WGN tracking error allows for tighter protection levels compared with colored noise.The large SS protection levels at small values of tracking error σ t are due to correspondingly large values of N min , which, in turn, lead to large coasting errors σ C k (with k N min = ).For colored noise, a shorter time constant provides a tighter protection level.Moreover, it can be seen that for a decimeter-level tracking error magnitude, the maximum protection level converges regardless of whether the tracking error is time-correlated.
When implemented together, the CPI and SS monitors address all of the limitations of the original CI monitor reported by Tanil et al. (2017), Tanil, Khanafseh, et al. (2018), Tanil, Jimenez, et al. (2018), and Kujur et al. (2019).The CPI monitor performs better than the CI monitor because its test statistic is optimized to

CONCLUSIONS
In this work, we introduced an optimal CPI monitor to detect spoofing by accumulating tracking error embedded in the spoofer's signal.We derived relationships between missed detection probability, tracking error magnitude, and monitor run time.We showed that even with decimeter-level tracking error, the monitor can detect spoofing with a low probability of missed detection in less than 1 min.We evaluated the performance of the CPI monitor for both white and time-correlated (colored) tracking error.To compute protection levels and detect short-duration spoofing, we proposed a complementary SS monitor to be implemented in sequential, overlapping windows to compare the integrated INS/GNSS position solution against an INS-only coasting solution.The INS-only coasting element also provides the capability to maintain positioning continuity after detection, albeit at lower accuracy, as the INS drifts.

APPENDIX A TIGHTLY COUPLED INS/GNSS ARCHITECTURE
Using IMU measurements, an INS provides a navigation state vector, which includes the aircraft position vector r with components x, y, z, velocity vector v with components u, v, w, and attitude φ, θ , ψ (Euler angles): An IMU consists of tri-axis accelerometers and gyroscopes to provide measurements of acceleration and body angular rate.In an INS, the IMU acceleration measurements are integrated once to obtain the velocity and then integrated again to obtain the position, whereas attitude is obtained by integrating angular rate measurements.These measurements have errors (biases and noise), causing the state estimate error to drift over time.In a tightly coupled INS/GNSS architecture, a KF uses raw GNSS code and carrier phase measurements to estimate and correct the error in the drifting INS states to provide an integrated navigation solution.
An individual (scalar) IMU measurement  u has errors such as time-dependent biases and noise.Therefore, it is modeled as a "true" measurement u * , corrupted with a constant bias b c , a slowly varying time-dependent bias-like component b, and additive WGN η u , as represented in Equation (A2).The constant bias is usually specified as bias repeatability, and the additive WGN η u is commonly derived from specifications on the VRW for an accelerometer and angular random walk for a gyroscope: The time-dependent component of the bias b is modeled as a first-order GMRP with time constant τ b and driving WGN ν b .The driving WGN is derived from the IMU "bias instability" specifications: The bias dynamics are included in the process model by augmentation of bias states x bias to the aircraft states.Thus, for three different IMU axes, the bias states for both acceleration and angular rate measurements are shown in Equation (A4).Equations ( A1) and (A4) show all of the nominal states that are propagated to obtain the INS navigation solution: We assume that an en route aircraft utilizes single-frequency GPS measurements without any differential corrections; however, this concept is also applicable to dual-frequency multi-constellation GNSS, terminal, and precision approach scenarios.Equation (A5) shows a simplified GPS measurement model in which the code measurement ρ for each satellite is composed of the true range r, satellite and receiver clock biases dt sv and dt rc , code ionospheric delay I ρ , code tropospheric delay T ρ , code multipath m ρ , and receiver code thermal WGN Q where c is the speed of light in vacuum and λ is the carrier wavelength.
All GPS errors must be included in the measurement in order to be utilized in the KF.The satellite clock offset cdt sv is corrected based on the clock parameters broadcast in the navigation message.After the satellite clock offset correction has been applied, there are still residual errors caused by satellite clock and orbit ephemeris parameter uncertainty.These residual errors r sv are modeled (Gallon et al., 2020) as a first-order GMRP with a time constant τ r sv of 5 h subject to driving WGN ν r sv with a standard deviation of 1.8 m.Equation (A6) represents the first-order GMRP model for satellite clock and ephemeris residual errors: where w r rc and w r rc  are WGN inputs to the clock offset and clock offset drift rate, respectively.The variances of these WGN inputs are obtained using typical Allan variance coefficients of temperature-compensated crystal oscillator timing standards.The white phase noise ( ) h 0 and frequency random walk noise ( ) h 2 coefficients used are 2 10 19 u and 2 10 20 u , respectively.For ionospheric delay, we use the ionospheric correction T iono from the Klobachaur model, which results in residual errors r i , as modeled by SC-159 (2009), with a standard deviation given by Equation (A8): where F pp is the obliquity factor and τ vert is calculated given the geomagnetic latitude (SC-159, 2009).Because the ionospheric delay is a slowly changing error, it is modeled as a first-order GMRP with a time constant of 40 h (Appendix B) and driving WGN ν r i as follows: The troposheric delay is corrected with the correction model specified in SC-159 (2009), and the residual errors r t in the zenith direction are modeled as a first-order GMRP with a time constant of 20 h and a standard deviation of 0.09 m (Gallon et al., 2021).The zenith error is converted to a slant ranging error using the elevation-dependent mapping function reported in SC-159 (2009).Equation (A10) shows the first-order GMRP model of the tropospheric residual error r t : where ν r t is the driving WGN for zenith tropospheric residual errors.Being time-correlated, the multipath is modeled as a first-order GMRP with a time constant τ m of 25 s and driving WGN ν m (SC-159, 2009): The standard deviation for code multipath error is 5 m (SC-159, 2009), and for carrier multipath error.we assume the standard deviation to be 0.02 m (Misra & Enge, 2012).The receiver code thermal noise standard deviation is taken to be 0.36 m (SC-159, 2009), and the carrier thermal noise standard deviation is assumed to be 3 mm (Misra & Enge, 2012).
Constant carrier phase cycle integer ambiguities, along with all of the above-mentioned residual error states, are included in the GNSS measurement error states: where n is the number of satellites in view.The final state vector of the INS/GNSS system is as follows:

B IONOSPHERIC RESIDUAL ERROR MODEL
Figure B1 shows a histogram of the ionospheric rate based on 24 h of data from Walter et al. (2004).The data are boundable by a Gaussian distribution with a mean of 3.2 10 2 u mm/s and a standard deviation of 0.96 mm/s.We model the

FIGURE B1
Histogram of rates of total electron content change for a quiet solar maximum day (Walter et al., 2004) ionospheric residuals using a first-order GMRP such that the obtained rate matches this Gaussian distribution.After utilizing Equation (A8) from SC-159 (2009) to obtain the standard deviation for the first-order GMRP, it was determined that the time constant that produces the best approximation of the ionospheric rate distribution in Figure B1 is approximately 40 h.

C EFFECT OF TRACKING ERROR ON COVARIANCE OF INNOVATIONS
We define the spoofer's uncertainty in the target's position as the tracking error ν t and use the superscript s for terms that are otherwise related to the spoofer or spoofed signal.At time k, the spoofer estimates the user position with tracking error, given as follows: where x y and z k k k , , give the true user position.We conservatively assume that there is no latency in the spoofer's tracking and counterfeit signal generation and neglect additive errors due to thermal noise and multipath in the counterfeit signal transmission.We also concede to the spoofer the ability to replicate the true signal-in-space errors seen by the target-namely, satellite orbit and clock errors and atmospheric effects.We also ignore the contribution of the deliberate deviations added by the spoofer to the counterfeit signal.Relaxing any of these conservative assumptions and concessions would only make the spoofing easier to detect at the target and the analysis less complicated.We ignore these effects in what follows and concern ourselves only with the impacts of tracking error in a single, arbitrary spatial direction.
The spoofer projects the counterfeit state vector x k us into the range domain:   (0, ) 2 is the tracking error along the spatial direction corresponding to the unit vector u in the state space.This vector can be re-expressed in terms of the "true" range vector, z k , as follows: The contribution of the tracking error to the innovations is then given by the following: Rearranging, we have the following relation: The covariance of the innovation vector with the presence of tracking error can be written as follows:

D EFFECT OF TRACKING ERROR ON TIME CORRELATION OF INNOVATIONS
The derivation of the Neyman-Pearson optimal test statistic relies on the assumption that the random variables forming the test statistic-the position-domain innovation in our case-are time-independent.To verify the independence of innovations over time, we generate sample auto-correlation functions (ACFs) and estimate the associated time constants.In the absence of tracking error, independence is guaranteed as long as the nominal error models in the KF are accurate.However, we also need to verify whether the assumption is valid if WGN tracking error is present.Figure D2 shows a sample ACF of carrier phase innovations for the case of WGN tracking error with V t 10 cm using a navigation-grade IMU.The time constant (measured at the 1/e point) of the sample ACF was 0.21 s, which is less than the sample interval of 0.5 s, thus confirming that the innovations are white.
We expect that a lower-quality IMU would have a greater contribution toward the time correlation of innovations because the Kalman gain would give a higher weight to the GNSS measurements than in the high-quality IMU case.Consequently, the tracking error will affect the post-measurement states x to a greater degree as well as the subsequent "predicted" measurement Hx used to generate the innovations.To confirm that the innovations are still white, even for lower-grade IMUs, we evaluate sample ACFs of the carrier phase innovations for navigation-, tactical-, and automotive-grade IMUs. Figure D2 shows the results for the case of WGN tracking error with V t 10 cm.It is clear from the figure that even for an automotive-grade IMU, the innovations still remain white, with a time constant of 0.24 s.

E COVARIANCE OF SS TEST STATISTIC
Let us consider an SS monitor that starts at time epoch 1.The KF time and measurement update equations are as follows: where, for any time epoch k, x k and ˆk x represent the states after the KF time and measurement updates, respectively.Φ k is the state transition matrix, L k is the Kalman gain, z k is the GNSS measurement, and H k is the observation matrix.The associated estimation equations can be determined as follows: where, for any time epoch k, e k and ˆk e represent the state estimate error after the KF time and measurement updates, respectively.I is the identity matrix, w k is the process noise vector, and v k is the measurement noise vector.For the SS monitor, parallel INS-only coasting is initiated with the state propagation equation: We can write the propagation equation for INS-only coasting as follows: The associated KF error propagation equations are as follows: The above equation can be further expanded: As before, this expression can be further expanded:

FIGURE 1
FIGURE 1 Performance of the CPI monitor versus the CI monitor for V t 2 cm (left) and for two different tracking error magnitudes (right)

FIGURE 2
FIGURE 2 CPI probability of missed detection P MD versus tracking error σ t and monitor run time N

FIGURE 6
FIGURE 6 Probability of missed detection P MD versus monitor run time N for WGN and colored noise tracking error The sample interval is 2 Hz; thus, N = 1 corresponds to 0.5 s.

FIGURE 8
FIGURE 8 Probability of missed detection ( ) P MD versus monitor run time (N) for different satellite geometries

FIGURE 9
FIGURE 9 Illustration of SS sequential monitor windows with increasing protection level (PL)

FIGURE 11
FIGURE 11 Maximum protection level (PL) versus tracking error magnitude for an example satellite geometry

U
th ( ).Similarly, the carrier phase measurement OI for each satellite is composed of the true range r, satellite and receiver clock bias dt sv and dt rc , carrier ionospheric delay I φ , carrier tropospheric delay T φ , carrier phase multipath m φ , carrier phase cycle integer ambiguity N φ , and receiver carrier thermal WGN Q I th .The code ionospheric delay I ρ is of the same magnitude as the carrier ionospheric delay I φ , and the code tropospheric delay T ρ is of the same magnitude as the carrier tropospheric delay T φ : offset cdt rc is compensated by a constant clock offset drift rate model.The clock offset state r rc is modeled to drift with a constant rate  r rc over time, as shown by Equation (A7):

FIGURE D2 Normalized
FIGURE D2Normalized ACF plots of the carrier phase innovations for different IMU grades with WGN tracking error (V t 10 cm) any time epoch k, x c k represents the state after the time update in INS-only coasting.Here, we also assume that the state transition matrix does not differ from that of the KF.The time update error propagation for the first time epoch of INS-only coasting is given as follows: any time epoch k, e c k represents the state estimate error after the time update during INS-only coasting.Similarly, for the second time epoch, we can write the KF time and measurement update equations as follows: The test statistic for the SS monitor at time k is defined as follows: E22) and (E25), we know that e 1 and e c 1 are the same for the first time epoch; hence, we have the following: the KF equations, we have the following expression for any time epoch k: . Because the SS monitor observes the instantaneous error between the tightly coupled INS/GNSS KF position solution and INS-only coasting solution and because the coasting covariance P c drifts over time, the detection threshold must