Impact Assessment of GNSS Spoofing Attacks on INS/GNSS Integrated Navigation System

In the face of emerging Global Navigation Satellite System (GNSS) spoofing attacks, there is a need to give a comprehensive analysis on how the inertial navigation system (INS)/GNSS integrated navigation system responds to different kinds of spoofing attacks. A better understanding of the integrated navigation system’s behavior with spoofed GNSS measurements gives us valuable clues to develop effective spoofing defenses. This paper focuses on an impact assessment of GNSS spoofing attacks on the integrated navigation system Kalman filter’s error covariance, innovation sequence and inertial sensor bias estimation. A simple and straightforward measurement-level trajectory spoofing simulation framework is presented, serving as the basis for an impact assessment of both unsynchronized and synchronized spoofing attacks. Recommendations are given for spoofing detection and mitigation based on our findings in the impact assessment process.


Introduction
Global Navigation Satellite System (GNSS) spoofing is a technique to trick the victim receiver into generating an erroneous position fix and/or clock offset, by deliberately broadcasting legitimate-appearing false satellite signals [1,2]. Civil GNSS service is vulnerable to spoofing due to the open structure and low power of GNSS satellite signals [1]. GNSS spoofing is difficult to detect and may result in more serious situations than jamming, since the user may not be aware of it [3].
GNSS spoofing is not a new topic. When Global Positioning System (GPS) was first declared full operational in 1995, researchers from MITRE examined GPS spoofing and its countermeasures [4]. The Volpe report, released by US Department of Transportation clearly pointed out the potential for spoofing attacks [1]. However, the GNSS community gave little attention to this threat, until a portable GPS spoofer was successfully developed and demonstrated [5]. Several proof-of-concept spoofing tests have been carried out against unmanned aerial vehicles [6,7], super yachts [8], cars and smartphones [6,9]. Besides these filed demonstrations, the Iran-US RQ-170 incident and the recently reported GPS problem in central Moscow [10] and the Black Sea [11], have further intensified the interest in this area. With open source GPS signal simulators available online and the fast developing software-defined radio technology, GNSS spoofing has become "not only feasible but also affordable" [12].
In response to the emerging GNSS spoofing threats, many contributions have been made, which can be divided into spoofing implementation, assessment, detection and mitigation. (1) Spoofing implementation focuses on the mechanisms of spoofing attacks, serving as the basis for impact analysis and validation of anti-spoofing methods. Several research groups have implemented the so-called receiver-spoofer as defined in [5], while the others rely on simulators/repeaters for spoofing simulation.
(2) Spoofing impact assessment evaluates the spoofing effects on GNSS receivers and systems that depend on, or relate to, GNSS. (3) Research on spoofing detection and mitigation, which take up the majority of the • Based on the authentic and spoofed signal model, a comparison between unsynchronized and synchronized GNSS spoofing attacks is given. A framework for a measurement-level trajectory spoofing simulation is proposed, which simplifies the process of impact assessment for integrated navigation systems.

•
We systematically analyze the impact of spoofing on the integrated navigation filter's error covariance, innovation and inertial sensor bias estimation, revealing how the conventionally used Kalman filer responds to spoofing attacks. • According to the impact assessment, we make recommendations for the cautious use of (1) error covariance for integrity monitoring, and (2) calibrated inertial sensors for pure INS solutions. Spoofing detection methods based on innovations and inertial sensor bias monitoring are suggested.
The remainder of this paper is organized as follows. Section 2 briefly introduces the INS/GNSS model used in this paper; Section 3 introduces basic GNSS spoofing attack modes and presents two methods for measurement-level spoofing simulation; Section 4 focuses on the impact assessment analysis and our recommendations for spoofing detection and mitigation; Section 5 verifies our findings through simulations, and the conclusions are given in Section 6.

INS/GNSS Integrated Navigation Model
In this paper, for simplicity, the loosely coupled INS/GNSS integrated navigation system is considered. As we focused on the Kalman filter's error covariance matrix, innovation and inertial sensor bias estimation, the analytical derivation and analysis in the following sections are also applicable to the tightly coupled system. The underlying INS/GNSS integrated system is a nonlinear time-varying system, so the extended Kalman filter was used for the integration. The error state, instead of the total state, of the navigation system was chosen as the state vector, which is defined as [22] T represent the gyro and accelerometer bias vectors, respectively. Note that the subscripts E, N, U represent the east, north and up components in the navigation frame, respectively. The subscripts, R, F, U represent the right, forward and up components in the body frame, respectively. After the linearization, the system dynamic model and measurement model can be written as [22] .
where F(t) is the system matrix; w(t) is the process noise vector; z(t) represents the difference between the position solution of INS p INS and GNSS p GNSS ; H(t) is the measurement matrix; and v(t) represents the noise of GNSS solutions. The system matrix, F(t), takes the form [23] where F ij is a 3-by-3 matrix; C n b is a 3-by-3 body frame to navigation frame rotation matrix; and 0 3 is a 3-by-3 zero matrix. The details of F ij are given in Appendix A [23]. The measurement matrix is a constant matrix and is defined as H(t) = 0 3 0 3 I 3 0 3 0 3 , in which I 3 is a 3-by-3 identity matrix.

GNSS Spoofing Attacks
A brief illustration of a GNSS spoofing attack is given in Figure 1. Assuming that the spoofing attack starts at t Init , the received GNSS signal, y(t), at the receiver can be expressed as [2], where y a (t) and y s (t) are authentic and spoofed signals, respectively. n(t) is the receiver noise. Taking GPS L1 coarse/acquisition code spoofing as an example, the received authentic signal y a (t) is given by The received spoofed signal y s (t) has a similar form: where the superscript (i) represents the ith satellite, the subscripts a and s represent the authentic and spoofed signals, respectively. N is the number of visible satellites; P is the signal power; D is the navigation data bit stream; C is the pseudo random noise (PRN) code sequence; τ is the signal propagation delay; f 1 is the GPS L1 signal frequency; f d is the Doppler shift; and θ 1 is the received initial carrier phase. There are several essential features of GNSS spoofing attacks [2,24]. The spoofing signals must have exactly the same PRN code sequence and signal frequency to those of the authentic signals. The number of spoofed satellites is generally equal to the number of authentic signals; if not, the introduction of inconsistency may trigger the receiver autonomous integrity monitoring algorithm. The structure of the navigation data bit stream should be the same as that of the authentic stream, but the content may be manipulated to introduce the intended deception. The received initial carrier phase is generally unknown and it is practically impossible for the spoofer to align with the authentic, unless a very precise (cm-level) relative position between the spoofer and the target receiver is known. The Doppler shift of the spoofed signal is not necessarily the same as the authentic signal; however, it must be consistent with the spoofer's own code phase variation.
GNSS spoofing attacks can be divided into unsynchronized and synchronized spoofing, based on whether the spoofed signals are time synchronized (code-phase aligned) with the authentic ones [16]. This classification highlights the error characteristics of spoofing attacks on raw GNSS measurements, which has an analogy to the step and ramp faults in integrity analysis. If the spoofer cannot know the target receiver's position accurately (meter-level to maximum of a half code chip) [16], only unsynchronized spoofing attacks can be carried out. For synchronized spoofing attacks, the spoofer accurately knows the target receiver's real-time position, which makes it possible for the spoofed signals' delay and the Doppler shift to be consistent with the authentic signal at the target receiver end [25]. A comparison between these two spoofing attack modes is given in Appendix B. Understanding the characteristics of different spoofing attacks is essential to achieve high fidelity spoofing simulation for impact assessment and defense verification.

Measurement-Level Spoofing Attack Simulation
To analyze the signal-level response of the GNSS receivers to spoofing attacks, a sophisticated signal-level spoofing simulator is necessary. However, for an impact assessment on INS/GNSS integrated navigation systems focusing on the Kalman filter performance with spoofed GNSS measurements, a measurement-level simulator is sufficient. The measurement-level simulation of spoofing attacks has the advantage of easy and rapid implementation, avoiding the complexity and high cost of setting up live-sky GNSS spoofing tests which must be authorized. The measurement-level simulator generates raw GNSS observations based on given a constellation, error model and predefined trajectory profile (each epoch with a seven-dimensional noise free PVT solution). Using the synchronized spoofing attack as an example, we show how the pseudorange measurements of the target receiver are affected, which also shows how the spoofer or spoofing simulator construct pseudorange measurements for a synchronized spoofing attack. While for unsynchronized spoofing, the spoofing simulator can generate arbitrary measurements without considering the relationship with the authentic signals.
A raw pseudorange measurement for the ith satellite ρ (i) at time t is modeled as where τ (i) is the signal transmission delay; c is the speed of light; and δt u and δt (i) are the receiver and satellite clock offset, respectively. For authentic signals, the signal delay, τ (i) a , consists of the delay caused by the geometric range, r Considering the receiver noise, n (i) ρ,a , the authentic pseudorange is obtained with For spoofed signals, the signal delay, τ where r (i) s is the geometric range from the spoofer to the satellite, and r s→u is the geometric range from the spoofer to the target receiver (common for all satellites). ∇τ (i) proc and ∇τ (i) ctrl are the signal processing delay and controlled signal delay, respectively. Assuming a common atmospheric effect for the spoofer and target receiver, Equation (10) can be written as where ∇τ (i) s is the additional signal delay introduced by the spoofer at the target receiver end. For a receiver/spoofer with the ability to predict the navigation bits, δt (i) adv_ctrl is introduced to represent the advanced prediction time by the spoofer to compensate for additional signal delays. Then, we get Assuming that ∇τ (i) proc is the same for all satellites, the superscript of ∇τ (i) proc will thus be omitted thereafter. In ∇τ (i) s , both the transmission delay, r s→u /c, and the signal processing time, ∇τ proc , will be estimated as the clock offset. Finally, we obtain where δt u,s and τ s (i) are the additional clock offset and signal delay introduced by the spoofing attacks.
While ∇τ proc can be calibrated by the spoofer [26], r s→u , r For measurement-level spoofing simulation, we define the authentic and spoofed trajectory profiles as tr a (t) and tr s (t), respectively. The two trajectories satisfy where s(t) is the desired spoofing profile which can be further defined as s(t) = α(t) + b. The term α(t) can be any time-related function and b is a constant vector. With α(t Init ) = 0 and b = 0, synchronized spoofing attacks can be simulated. As a good start, α(t) can be a simple ramp function similar to integrity analysis in the GNSS community. For unsynchronized spoofing, the term b is used to represent the jumps introduced to the authentic GNSS solutions. There are two ways to achieve the measurement level spoofing attack simulation: • When s(t) is determined, the components δt u,s and τ s (i) in Equation (13) can be calculated accordingly. They are added to the authentic measurements generated based on tr a (t) to construct the spoofed measurements. In this way, the critical parameters of the spoofer can also be simulated.

•
If we only focus on the INS/GNSS integrated navigation system, a simpler method can be used without direct calculation of δt u,s and τ s (i) . This is done by directly feeding tr a (t) and tr s (t) to the GNSS measurement-level simulator as two independent trajectories. A switch from tr a (t) to tr s (t) during the simulation can easily introduce the spoofed measurements to the integrated navigation system.

Error Covariance
The diagonal elements of the Kalman filter error covariance matrix, P k , represent the error variance of each state estimation when the filter is properly modeled. They are often used to evaluate the performance of the integrated navigation system. In face of GNSS spoofing attacks, there is a necessity to analyze its impacts on the Kalman filter P k calculation.
A standard Kalman filter implementation is given in Figure 2 with spoofing attacks introduced. GNSS spoofing attacks are injected directly into the state filtering loop by adding a spoofing profile to z k . For a typical spoofing attack, a positioning error of several tens of meters to several kilometers is sufficient to cause serious problem for safety critical applications. It should be noted that, in order to analyze the effects solely caused by the spoofing attacks, the system dynamic/measurement model (including the associated model parameters such as the process noise covariance matrix, Q k and the measurement noise covariance matrix, R k ) and the filter initialization process, are assumed to be exactly the same under spoofed and authentic conditions in the impact assessment analysis. This ensures that all the other factors that may influence the Kalman filter's performance are excluded. As shown in Figure 2, the only link between the state filtering loop (left) and the gain calculation loop (right) is the Kalman filter gain matrix, K k . Intuitively, this is a one-way connection which means that the left loop is affected by the right but not vice versa. For linearized Kalman filter implementation, in which the system equation is linearized along the predefined reference trajectory, the gain calculation loop (which can be done off-line) is completely independent to the state filtering loop. For the more commonly used extended Kalman filter, there is a minor difference between the open-loop feedforward and closed-loop feedback mechanisms. If the feedforward integration is used, in which only the outputs of the integrated navigation system are corrected, the pure INS-based reference trajectory remains unaffected, and thus, the state filtering loop is also completely independent of the gain calculation loop. However, for the feedback integration, the error states estimated by the Kalman filter are used to correct the velocity, position, attitude, and inertial sensor measurements of the inertial navigation system, which means that the reference trajectory is corrected and updated periodically, and thus, the state filtering loop will affect the gain calculation loop through the calculation of Φ k/k−1 . However, the changes in position induced by typical spoofing attacks do not significantly influence the Φ k/k−1 calculation. A simple and intuitional numerical demonstration is given in Appendix C, in which an assumption of a 5 km position error is introduced to each of the three dimensions. The maximum relative change of each element of Φ k/k−1 is less than 2.5%. This is because a position change of several kilometers is only a small perturbation compared to the earth radius for the Φ k/k−1 calculation, as shown in Appendix A. Generally speaking, the gain calculation loop is weakly correlated with the state filtering loop.
As explained above, the spoofing profile introduced to the state filtering loop will have little impact on the Φ k/k−1 calculation, so P k will basically remain unchanged under the assumption that the other terms (Q k , R k , H k ) and the filter initialization are exactly the same under spoofed and authentic conditions. This can be verified from the simulations in Section 5. In other words, P k can no longer reflect the filter's performance under spoofing attacks. The basic premise of using P k for performance evaluation is proper modeling of both system dynamics and observations. Under spoofing attacks, the measurement model is considered to be incorrect, because the spoofing profile is not modeled and the constant R k cannot reflect the error of the spoofed GNSS measurements. The measurement model used under authentic conditions is no longer valid when the spoofing attacks are introduced. This is the main reason why P k is unreliable and cannot reflect the filter's performance under spoofing attacks. Therefore, P k should be used cautiously to evaluate the navigation performance, taking spoofing attacks into account.
In civil aviation, the commonly used integrity monitoring algorithms, Honeywell Inertial GPS Hybrid (HIGH) and Autonomous Integrity Monitored Extrapolation (AIME), calculate the horizontal protection level based on P k [27]. Under GNSS spoofing attacks, if the spoofed signals are not detected and removed in time, the horizontal protection level derived from P k is unreliable, which may cause potential integrity risk. Here, we recommend that the integrity algorithm designers reconsider and evaluate the integrity monitoring methods that rely on or relate to P k , with spoofing attacks taken into consideration.

Kalman Filter Innovation
The Kalman filter innovations are often used to construct the fault detection statistics for INS/GNSS integrated navigation systems. It is necessary to analyze the effects of spoofing attacks on the Kalman filter innovations to gain a better understanding about spoofing detection methods that are based on them. The discrete system dynamics and measurement model can be written as [22].
The Kalman filter time and measurement update process are implemented in accordance with Figure 2. The Kalman filter innovation, υ k , is defined as [28] Under normal conditions, when the filter approaches a steady state, the innovation has zero expectation and known covariance [28], that is The test statistic of a Chi-squared test for fault detection in the INS/GNSS integrated navigation system at epoch k has the simplest form as [22] If there is no failure, the test statistic obeys a central Chi-squared distribution with m degrees, where m is the dimension of the measurement vector. The detection threshold is determined given a constant false alarm rate based on the Chi-squared distribution. The Kalman filter innovations can be used to construct a variety of test statistics, and Equation (20) gives the simplest and most commonly used one. The mean of the innovations can be tested and the innovations within a time window can be accumulated or averaged to build different forms of Chi-squared tests [19,23,29].
GNSS spoofing detection is a complex fault detection problem, in which the spoofing attacks can be modeled as drifts or abrupt changes in the sensor measurements. Change in the expectation of the innovation sequence under a drift fault has been investigated in the pressure sensor drift detection problem of a nuclear power plant [30]. Here, we generally follow the derivation given in [30] but extend the fault profile from a simple ramp manner with a fixed rate to a general spoofing attack vector to clearly show the statistical property of the innovation sequence under spoofing attacks. Assume that the spoofing vector, S n (n = 0, 1, 2, 3 . . . ), is introduced from t Init . For simplicity, we use the new subscript, n, to represent the time epoch in a spoofing attack. Now the measurement model becomes The Kalman filter innovation under spoofing attacks turns to Define the one-step prediction error at time n as Then, the expectation of the innovation under spoofing attacks is The measurement update equation can be rewritten aŝ x n =x n/n−1 + K n H n (x n −x n/n−1 ) + K n (V n + S n ).
Based on the system dynamic model, We obtain the one-step prediction error at time n + 1 by subtracting Equation (28) from (29): Then, the expectation of x n+1 is Given a fault free initial condition, E x 0 = 0, the expectation of x n can be obtained from where Finally, we obtain the analytical expectation of the Kalman filter innovation under spoofing attacks: From Equation (34), it is obvious that the spoofing profile at the current epoch has a direct effect on the expectation of the current innovation, while the spoofing profile history has indirect, but accumulated, impacts. If the spoofing profile of each epoch is known, Equation (34) can be solved recursively. As the expectation of the innovation changes under spoofing attacks, carrying out a statistical test directly on the mean of the innovation is a simple and straightforward spoofing detection method. For the error covariance of the innovation, if the spoofing vector does not include additional measurement noise, the error covariance will remain the same as that defined in Equation (19), that is P υ n,s = H n P n/n−1 H T n + R n .
Meanwhile, the test statistic in Equation (20) does not obey a central Chi-squared distribution anymore, but changes to a non-central Chi-squared distribution with the non-centrality parameter calculated with [19,21] q n = E(υ n,s ) T P −1 υ n,s E(υ n,s ).
The change in the statistical property of the innovations under spoofing attacks leads to a series of innovation-based spoofing detection methods. Furthermore, Tanil [19,21] established a novel, worst-case spoofing profile aimed at maximizing the position error without being detected, based on a similar expression of Equation (34) with slightly different derivations. The worst-case spoofing profile was used instead of the typical ramp and step fault profile that is commonly seen in the GNSS community to evaluate the spoofing detection performance. A linearized Kalman filter implementation was used in [19,21], in which a predefined and known nominal trajectory greatly simplified the spoofing detection and evaluation problem. However, the application of the linearized Kalman filter for INS/GNSS integration is limited to very certain type of applications, like the precision approach, orbiting satellites or interplanetary travel. For the more commonly used extended Kalman filter integration without a predefined reference trajectory, with linearization carried out along the real-time estimation, the spoofing profile is tightly coupled with the Kalman filter real-time estimation, so the worst-case spoofing profile derivation is not applicable. For the general INS/GNSS integration navigation system, typical ramp and step fault profiles (corresponding to synchronized and unsynchronized spoofing attacks) are used in this paper for the simulation analysis. The existence and derivation of a more advanced or, so-called, worst-case spoofing profile that maximizes the integrity risk still needs further effort. The analytical expression derived here could serve as a basis for further investigations.

Inertial Sensor Bias Estimation
In addition to introducing errors in the position and velocity solutions, the spoofing attacks also affect the estimation of inertial sensor biases. Following the derivation in Section 4.2, as we have already obtained the expectation of the one-step prediction error, x n , and the expectation of the innovation, υ n,s , we can further obtain the expectation of the inertial sensor bias estimation error under spoofing attacks using the following derivation.
First, define the state estimation error at time n as The relationship of the state estimation error, x n , and the one-step prediction error, x n , can be obtained by combining Equation (37) and Equation (23). Then, we get The state estimation is given bŷ x n =x n/n−1 + K n (z n,s − H nxn/n−1 ) =x n/n−1 + K n υ n,s .
Substituting Equation (39) into Equation (38), we get As the expectations of x n and υ n,s have already been given in Equations (32) and (34), respectively, the expectation of the state estimation error can be obtained with where the term, Λ n,i , has been defined in Equation (33). Recall that the definition of the state vector is ; the last six elements of E( x n ) give the expectation of the error of inertial sensor bias estimation. As with the explanation given under Equation (34), the inertial sensor bias estimation is affected by the spoofing profile in an accumulated way. If the spoofing profile of each epoch is exactly known, Equation (41) can be solved. However, as we pointed out at the end of Section 4.2, when the extended Kalman filter is used, the spoofing vector will be tightly coupled with the real-time state estimation, which makes it difficult to build a quantitative relationship between a given spoofing profile and its impacts on the inertial sensor bias estimation without a full simulation analysis. We use a simple analytical demonstration below to give a qualitative analysis here. A simulation is given in Section 5 for more intuitional and quantitative demonstration.
As a common way to demonstrate and analyze the error propagation, we consider a simple static situation of the north channel for qualitative analysis of spoofing impacts on the inertial sensor bias estimation. The north position, north velocity and east misalignment angle error equation can be simplified with Equations (42)-(44), respectively [31]. .
where δr N , δv N and φ E are the north position error, north velocity error and east misalignment angle, respectively. ∇ N and ε E represent the north accelerometer bias and the east gyro bias, respectively. g and R represent the gravity and the earth radius, respectively. Assuming that the spoofer introduces a north position or velocity error, as shown in Equation (43), both the east misalignment angle and north accelerometer bias will be affected. Furthermore, when the east misalignment angle is influenced, the east gyro bias is affected through Equation (44). Generally speaking, a spoofing attack on the north will introduce estimation errors on the north accelerometer bias, east gyro bias, and east misalignment angle. Note that the above analysis only serves as a simple demonstration. The impacts will be much more complicated due to cross coupling of different errors in dynamic situations.
As the gyro and accelerometer bias estimations are affected by the spoofing attacks, spoofing attacks may be detected by monitoring the bias estimation of the inertial sensors from the Kalman filter. This can be done by directly setting upper and lower bounds on the estimated biases. Based on the authors' experience, the threshold should be set at least three to five times larger than a nominal value provided by the manufacturer. This method is generally conservative, and the user should take the risk that the inertial sensors suffer from performance degradation for other reasons rather than spoofing attacks.
After successful detection of the spoofing attacks, the integrated navigation filter should discard the GNSS solutions if spoofed signals are not removed in time. The navigation filter will output a pure INS solution instead. Conventionally, the estimated inertial sensor biases are used to compensate the raw gyro and accelerometer measurements in the pure INS solution when GNSS measurements are not available for jamming or blockage. In the spoofing scenario, however, the estimated inertial sensor biases are no longer reliable and should not be used. Meanwhile, the spoofing attacks can also affect the misalignment angle estimation, as illustrated above, which also leads to large, pure INS errors. To mitigate the spoofing impacts, the backtracking mechanism, which records the raw IMU measurements in a time window and starts the pure INS solution from a backward time, is recommended after successful spoofing detection. This mechanism has been applied to the initial alignment [32] and can be transferred to spoofing mitigation with minor modifications.

Simulation Analysis
To verify our analyses, we carried out a series of simulations based on the proposed measurement-level trajectory spoofing simulator presented in Section 3.2. We directly fed the spoofing and authentic trajectories to the GNSS measurement-level simulator and switched from authentic to spoofing trajectories during the simulation to introduce the designed spoofing attacks. A GNSS positioning engine was implemented separately into the simulation, with spoofed pseudoranges as inputs and position solutions as outputs for the loosely coupled integration. A flight trajectory lasting 1100 s near Xi'an, China was simulated as the authentic trajectory. The simulated position profiles are given in Figure 3. The parameters of the Kalman filter used in the simulation are listed in Appendix D. We introduced the spoofing attacks (in the longitude channel in all cases) in the straight-and-level flight period, starting at 800 s; the total period of spoofing was 300 s. For the GNSS simulation, we used the baseline 24-slot GPS constellation at an epoch of 00:00:00 on 1 July 1993 [33], which is also recommended in the RTCA/DO-229 MOPS for INS/GPS system evaluation. The aviation-grade INS and single-frequency GPS receiver were adopted in the simulation with the specifications listed in Table 1. All the following simulations were carried out under the same dynamic scenario as described above.   We further compared each diagonal element of P k under spoofed and authentic conditions for all of the estimated states in the scenario shown in Figure 4b at t = 1100 s, when about 2 nautical miles longitude error had been introduced. The absolute change in the square root of each diagonal element of P k was defined as where the subscripts, s and a, represent the spoofed and authentic conditions, respectively. i = 1, 2, 3, . . . , 15 represents the index of the diagonal element of P k . The relative change was defined as Both the absolute change and relative change are given in Table 2, which clearly shows that the changes in P k -indicated state estimation errors were nearly negligible. This verifies that spoofing attacks injected into the state filtering loop have little impact on Kalman filter error covariance calculation. Using P k to evaluate the integrated navigation system's performance and for integrity monitoring is unreliable under spoofing attacks.
Note: absolute and relative represent the absolute change, ∆sqrtPk, as defined in Equation (45), and relative change,δsqrtPk, as defined in Equation (46), respectively. Figure 5a shows the three dimensional innovations under a 0.5 m/s spoofing attack. The abnormal innovations were found only in the longitude channel, because the spoofing attacks introduce only longitudinal errors. Separately, tests on the mean of each component of the innovation vector have the potential to detect and distinguish which component is affected. In Figure 5b, the longitude channel innovations are given for different spoofing profiles, which shows that it may be easier to detect spoofing attacks with larger fault ramps. Also, there is only a limited spoofing detection window [34], as the Kalman filter dynamically tunes itself to track the spoofing profiles, and the means of innovations will approach to zero again in a new steady state. In addition to testing the means of innovations, the other innovation-based tests, like the snapshot Chi-squared test [22,23] or averaged/summed innovation test [19,23,29], can also be used. The impacts of spoofing attacks on the inertial sensor bias and misalignment angle estimation are illustrated in Figure 6a. Taking the gyro bias estimation as an example, the Kalman filter gives an estimation of the y-axis gyro bias of over 0.3 • /h under a 1 m/s synchronized spoofing attack, far beyond the 0.01 • /h specification. This gives an obvious sign that abnormal GPS measurements exist. Assuming that the spoofing attack is detected at 825 s, and the GPS measurements in the integrated system are directly disconnected, the horizontal error of pure INS grows quickly to about 150 m, as shown in Figure 6b. When the backtracking mechanism given in [32] (with 50 s backtracking) was used, the unaffected inertial sensor bias and misalignment angle estimation contributed to a better pure INS solution performance with position errors about 90 m. This shows the potential benefits of the backtracking mechanism for spoofing mitigation.

Conclusions and Future Works
Understanding the behavior of INS/GNSS integrated navigation systems under spoofing attacks is crucial to allow the development of effective spoofing detection and mitigation methods. In this paper, we analyzed the spoofing impacts with a focus on the Kalman filter of the integrated navigation system. Through theoretical analysis and simulations, we came to the conclusion that (1) Kalman filter state error covariance is not affected by spoofing attacks and thus, should be reconsidered for performance evaluation and integrity monitoring; (2) the statistical properties of innovations have changed, which leads to several potential statistical tests for spoofing detection; (3) the estimated inertial sensor biases are no longer reliable under spoofing attacks, which makes bias estimation range check and the backtracking mechanism promising for spoofing detection and mitigation.
What we touched on in this paper is only the tip of the iceberg; this paper opens many research opportunities for the integrated navigation community, taking spoofing attacks into consideration. Simulations in this paper give an intuitional view of the possibility of innovation-based and inertial sensor bias monitoring methods for spoofing detection. Detailed investigations of the recommended spoofing defenses with a comparison and maybe a combination of these defenses deserve further efforts.

Acknowledgments:
We would like to express our great appreciation to the anonymous reviewers for their careful reading of our manuscript and their insightful comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The details of each component of the system matrix, F(t) can be written as where ( * )× denotes the skew-symmetric matrix of the corresponding vector. R Mh and R Nh are the meridian and transverse radii of the curvature plus the height.
 is the body frame specific force vector, consisting of right, forward and up components.
ω ie is the constant earth rotation rate (=7.292115 × 10 −5 rad/s). ω n ie is the earth rotation vector resolved into the local navigation frame, defined as ω n en is the rotation vector of the navigation frame with respect to the earth-fixed-earth-centered frame represented in the local navigation frame, defined as Appendix B

Appendix C
To give a simple and intuitional demonstration of the effect of introduced position error on the system state transition matrix Φ calculation, we take a snapshot of the simulation scenario presented in Section 5. Under the authentic condition, at t = 1100 s, the snapshot system state is given as Assuming that 5 km position errors (denoted by δL, δλ and δH) are directly added to the reference trajectory for all the three dimensions, elements of the state transition matrix, Φ i,j (element in the ith row and jth column), are calculated before and after the introduction of position errors. Let Φ i,j(a) and Φ i,j(s) represent the element under authentic and spoofed conditions, respectively. The elements of δΦ, representing the relative change of Φ, are defined and calculated as where δΦ i,j is the element in the ith row and jth column of δΦ. The first order solution, I + F(t)T s , is used to compute the transition matrix, where F(t) is the system matrix, as defined in Equation (3) and Appendix A, and T s is the propagation interval which equals 0.01 s in this demonstration. Here, take the calculation of Φ 8,4 as an example (corresponding to the term of sec L/R Nh in Equation (A7)). The relative change of Φ 8,4 is calculated as (A15) Note that only part of δΦ is given; all the omitted elements and the elements denoted by "-" are not influenced by introduced position errors (elements with relative change <1 × 10 −9 are also denoted by "-") .The superscripts, * and **, are used to indicate δΦ 8,4 for demonstration in Equation (A14) and δΦ 5,9 for the maximum relative change, respectively.

Appendix D
The parameters of the Kalman filter used in the simulation are listed here for reference: x 0 = zeros (15,1).
Note that Ts = 0.01 s is the propagation interval in the simulation. The parameters are set according to IMU and GPS specifications, as listed in Table 1. The units given above are just for easy interpretation. They should be changed to SI in practical implementations.