Tuning-Free, Low Memory Robust Estimator to Mitigate GPS Spoofing Attacks

IEEE CONTR OL SYSTEMS LETTERS, JUNE 2019, IN PRESS 1 T uning-Free, Lo w Memory Rob ust Estimator to Mitigate GPS Spoofing Attacks Junhwan Lee, Ahmad F . T aha, Nikolaos Gatsis, and David Akopian. Abstract —The operation of critical infrastructures such as the electrical power grid, cellphone to wers, and financial in- stitutions relies on precise timing provided by stationary GPS recei vers. These GPS devices are vulnerable to a type of spoofing called Time Synchronization Attack (TSA), whose objective is to maliciously alter the timing provided by the GPS receiver . The objective of this paper is to design a tuning-free, low memory rob ust estimator to mitigate such spoofing attacks. The contribution is that the proposed method dispenses with several limitations found in the existing state-of-the-art methods in the literature that r equire parameter tuning, availability of the statistical distributions of noise, real-time optimization, or heavy computations. Specifically , we (i) utilize an observer design for linear systems under unknown inputs, (ii) adjust it to include a state-correction algorithm, (iii) design a realistic experimental setup with real GPS data and sensible spoofing attacks, and (iv) showcase how the pr oposed tuning-fr ee, low memory rob ust estimator can combat TSAs. Numerical tests with real GPS data demonstrate that accurate time can be provided to the user under various attack conditions. Index T erms —Robust state estimation, observer design, GPS spoofing, time synchronization attacks, low memory estimation. I . I N T R O D U C TI O N A N D M OT I V AT I O N T HE Global Positioning System (GPS) is widely uti- lized in an abundance of applications. The study [ 1 ] in particular emphasizes ho w critical infrastructures such as communications, the power grid, transportation, and even financial services can be disrupted if the integrity of the GPS is compromised. Since most systems rely on non-encrypted civilian GPS signals [ 2 ], the GPS is vulnerable to intentional attacks. There are two types of deliberate attacks on GPS: jamming and spoofing [ 3 ]. While jamming completely blocks signal reception by transmitting high power noise, spoofing changes the transmitted signal or data to deceive the GPS receiv er . V arious experiments have sho wn that different types of spoof- ing attacks such as data lev el spoofing, signal le vel spoofing, delaying, and record-and-reply attacks [ 4 ]–[ 6 ] could af fect off- the-shelf GPS recei vers. A T ime Synchronization Attack (TSA) is a particular class of spoofing attacks on stationary GPS receiv ers that pro- vide precise timing in various applications, including phasor measurement units (PMUs), cellphone towers, and financial institutions [ 7 ], [ 8 ]. The objectiv e of the attack is to mislead the time estimated by the recei ver which, in the case of the The authors are all with the Department of Electrical and Computer Engineering at the Univ ersity of T exas, San Antonio. { junhwan.lee,ahmad.taha,nikolaos.gatsis,da vid.akopian } @utsa.edu. W e acknowledge the financial support by the National Science Foundation under Grant ECCS-1719043, and the suggestions made by the editor and re viewers. power grid for example, can disrupt the reliable monitoring of the grid’ s cyber -physical status. Countermeasures against spoofing attacks have been pro- posed in [ 9 ]–[ 11 ] and include techniques that rely on multiple GPS recei vers [ 12 ], [ 13 ] or check the magnitude of error in the GPS data [ 9 ]; see [ 5 ] for a revie w of anti-spoofing techniques. Another approach to mitigate and detect spoofing TSAs is through robust, dynamic state estimation routines that are designed to deal with outliers and malicious attacks. The robust state estimation literature is indeed rich with two main classes of methods. The first class is based on robust observers and L yapunov theory which often does not assume any statistical distribution for unknown inputs or noise [ 14 ], [ 15 ]. The second class is based on Kalman filter and its deriv a- tiv es, that often assume statistical distrib ution of noise [ 16 ]– [ 18 ]. Rele v ant to the rob ust estimation problem under TSAs, a nov el anti-spoofing particle filter is devised to find the receiver position ev en under spoofing interference [ 10 ]. In our recent work [ 11 ], we develop a real-time optimization method to detect and mitigate TSAs using weighted ` 1 minimization. The aforementioned methods all require either tuning of some parameters, real-time optimization, or av ailability of the statistical distribution of noise. The objecti ve of this paper is to design a robust state estimator that combats TSAs while being endowed with the following properties: (i) It is a tuning-free method that does not require any training; (ii) it has low memory requirement in the sense no heavy computations are needed in real-time; and (iii) the designed tuning-free, low-memory robust estimator can correctly reconstruct the actual physical state of the receiver using a realistic testbed. W e note here that our objective is not to develop a generalized theory for robust estimation, but rather to build on the recent theoretical advancement in this field and adapt it to the GPS spoofing problem through a realistic testbed. The rob ust state estimation method presented in this paper is an adaptation of the method in [ 19 ]. The motiv ation for using this estimator is pro vided in great detail in the next sections; the paper’ s organization is gi ven next. Section II presents the dynamic modeling of bias and drift in GPS receivers and sho wcases how TSA attacks can be modeled and designed to mislead standard state estimators. Section III presents a robust state estimator in addition to a state correction and TSA reconstruction routine. Section IV concludes the paper with realistic numerical tests on real GPS data. For reproducibility , the data used in the numerical tests and the results are all provided through a Github link. 2 IEEE CONTR OL SYSTEMS LETTERS, JUNE 2019, IN PRESS I I . D Y N A M I C M O D E L I N G U N D E R G P S S P O O FI N G The primary goal of GPS localization is to accurately estimate the position, velocity , clock bias, and clock drift of the receiver in ev ery time step—con ventionally referred to as the p osition, v elocity , and t ime (PVT) solution. The position of the GPS receiver (user) in Earth Centered Earth Fix ed coordinates is denoted by p u = [ x u , y u , z u ] > . T o estimate the receiv er’ s location and velocity , the GPS exploits the kno wn location of satellites and the distance from each satellite to the receiv er . Let N denote the num- ber of fixed satellites visible by the recei ver , then p n = [ x n ( t n ) , y n ( t n ) , z n ( t n )] > for n = 1 , 2 , . . . , N is the satellite position at the time of transmission t n . Also, we consider that t R models the arriv al time of transmitted signal at the receiv er . The approximate distance between each satellite and the receiv er can be written as ρ n = c ( t R − t n ) , where c is the speed of light and ρ n is commonly known as the pseudorange [ 20 ]. The pseudoranges differ from the true distances because t R and t n are offset by clock biases denoted by b u and b n . This relationship is established as t R = t GPS R + b u , t n = t GPS n + b n (1) where t GPS R defines the reception time in the absolute GPS time, while t GPS n is the signal transmission time in GPS time. Then, if d n represents the true distance from the receiv er to each satellite, it holds that d n = c ( t GPS R − t GPS n ) using the unbiased transmission and reception times. Alternatively , d n can be expressed by taking the 2-norm of position difference between the satellite and the receiv er , given by d n = || p n − p u || 2 . The following pseudorange equation is generated by combining the pre vious two equations for d n [ 20 ] ρ n = k p n − p u k 2 + c ( b u − b n ) +  ρ n (2) where  ρ n captures atmospheric ef fects and recei ver noise. In addition to the pseudoranges, the GPS receiv er can also measure the rate at which the pseudoranges vary over time, denoted by ˙ ρ n , and called pseudorange rate. The pseudorange rates are expressed in terms of the satellite velocities v n and the user (GPS receiv er) velocity v u as ˙ ρ n = ( v n − v u ) > p n − p u k p n − p u k + ˙ b u +  ˙ ρ n (3) where ˙ b u represents the GPS recei ver clock drift and  ˙ ρ n is the noise. In ( 2 ) and ( 3 ), the unknown PVT variables ()user position ( p u ), user velocity ( v u ), clock bias ( b u ), and clock drift ( ˙ b u )] are usually computed using nonlinear weighted least squares. The random walk model captures the dynamics relating variables in ( 2 ) and ( 3 ) for stationary applications [ 20 ]. The following is the stationary random walk model:       x u [ k + 1] y u [ k + 1] z u [ k + 1] b u [ k + 1] ˙ b u [ k + 1]       =   I 3 × 3 0 3 × 2 0 2 × 3 1 ∆ t 0 1         x u [ k ] y u [ k ] z u [ k ] b u [ k ] ˙ b u [ k ]       + w [ k ] (4) where k is the time index; ∆ t is the time resolution; and w is noise in the system. Generally , stand-alone recei vers like those present in PMUs, use the Extended Kalman Filter (EKF) to estimate the PVT solution [ 20 ]. Since the receiv er is stationary , the position ( p u ) can be treated as known constant while the recei ver v elocity ( v u ) is known to be zero. Thus, the only variables to be estimated are in fact the clock bias and drift, b u [ k ] and ˙ b u [ k ] . Based on ( 2 ), ( 3 ) and the dynamic model in ( 4 ), the fundamental plant model is constructed as follows using ρ [ k ] = [ ρ 1 [ k ] , ..., ρ N [ k ]] > and ˙ ρ [ k ] = [ ˙ ρ 1 [ k ] , ..., ˙ ρ N [ k ]] > :  cb u [ k + 1] c ˙ b u [ k + 1]  = A  cb u [ k ] c ˙ b u [ k ]  + w [ k ] (5)  ρ [ k ] ˙ ρ [ k ]  = C  cb u [ k ] c ˙ b u [ k ]  + c l [ k ] +  [ k ] (6) where A =  1 ∆ t 0 1  , C =  1 N × 1 0 N × 1 0 N × 1 1 N × 1  c l [ k ] =               k p 1 [ k ] − p u [ k ] k − cb 1 [ k ] . . . k p N [ k ] − p u [ k ] k − cb N [ k ] ( v 1 [ k ] − v u [ k ]) > . p 1 [ k ] − p u [ k ] k p 1 [ k ] − p u [ k ] k − c ˙ b 1 [ k ] . . . ( v N [ k ] − v u [ k ]) > . p N [ k ] − p u [ k ] k p N [ k ] − p u [ k ] k − c ˙ b N [ k ]               and w [ k ] and  [ k ] represent process/measurement noise; vec- tor c l [ k ] is based on the known satellite position, velocity and clock characteristics—a time-v arying, kno wn quantity . Equations ( 5 ) and ( 6 ) can be written as x [ k + 1] = Ax [ k ] + w [ k ] y [ k ] = C x [ k ] + c l [ k ] +  [ k ] . (7) The state space model ( 7 ), howe ver , does not model poten- tial spoofing attacks. While many different physical spoofing mechanisms are devised to deceive the victim receiv er [ 6 ], time synchronization attack (TSA) is applied on the stationary GPS receiv er . In practical sense, TSA alters the timestamp estimate by inserting the spoofing signal into the authentic pseudorange signals: ρ s [ k ] = ρ [ k ] + s ρ [ k ] , ˙ ρ s [ k ] = ˙ ρ [ k ] + s ˙ ρ [ k ] . where s ρ [ k ] and s ˙ ρ [ k ] denote the spoofing attacks, and ρ s [ k ] and ˙ ρ s [ k ] are the spoofed measurements. Specifically , there are two dif ferent types of TSAs according to the shape of s ρ . While T ype I attack injects an abrupt signal, e.g., s ρ [ k > α ] = 8000 m where α indicates the initial time of attack, T ype II attack modifies the clock bias in gradual manner manipulated by s ρ [ k ] = s ρ [ k − 1] + s ˙ ρ [ k ]∆ t ; see [ 21 ]. The actual effect of each type of attack on the clock state is thoroughly revie wed in [ 22 ]. As an example, in order for the spoofing signal to be considered as intentional attack on a PMU, it has to satisfy certain conditions. According to the IEEE C37.118 Standard, the attack has to result in 1% total variation error , which is equiv alent to 26.65 µ s clock bias error , or 7989 m of distance equiv alent bias error in order for the attack to be considered LEE ET AL. ; TUNING-FREE, LOW MEMOR Y R OBUST ESTIMA TOR TO MITIGA TE GPS SPOOFING A TT ACKS 3 infringing [ 23 ]. These types of spoofing attacks—regardless of their physical mechanism—impact the state dynamics as x [ k + 1] = Ax [ k ] + d [ k ] where d [ k ] = [ d 1 [ k ] d 2 [ k ]] > models and lumps TSAs and any process noise. Concrete examples of TSAs are given in Section III . In our pre vious work [ 11 , Section III], we show how specific forms of d [ k ] can mislead the receiver . I I I . R O BU S T S T AT E E S T I M A T O R In this section, we present a state estimation algorithm that is endowed with the following properties: (i) It is a tuning- free method that does not require any knowledge of noise distribution, initial parameters or states, or other coef ficients; (ii) it has low memory requirements in the sense no heavy computations are needed which is befitting to devices with limited computational power and limited internet connectivity; (iii) it is robust to GPS spoofing, time-synchronization attacks. A. GPS Clock Model and Estimator Dynamics The plant model under the spoofing attack based on the previous can be written as x [ k + 1] = Ax [ k ] + d [ k ] y [ k ] = C x [ k ] + c l [ k ] +  [ k ] (8) where x [ k ] ∈ R 2 represents the state v ector of clock bias and drift at time k ; y [ k ] ∈ R 2 N represents a single column vector of pseudoranges and pseudorange rates where N indicates the fixed number of visible satellites at e very time inde x; d [ k ] ∈ R 2 is the unknown spoofing attack applied to the bias and drift which also includes process noise; state-space matrices A , C , c l are discussed in the pre vious section. Consider no w a new modified state vector x m [ k ] ∈ R 2 which represents the state vector without spoofing attack and follows the following dynamics: x m [ k + 1] = x [ k + 1] − d [ k ] = Ax [ k ] . (9) The left-hand side of ( 9 ) essentially represents the original state vector considering that the spoofing attack d [ k ] is re- mov ed. The modified state vector x m [ k ] propagates through to y [ k ] . This yields: x m [ k + 1] = Ax m [ k ] + Ad [ k − 1] y [ k ] = C x m [ k ] + c l [ k ] + C d [ k − 1] . (10) The presented state estimator in this paper is an adaptation of the observer from [ 19 ] and follo ws the dif ference equation: R O B U S T E S T I M ATO R ˆ x m [ k + 1] = A ˆ x m [ k ] + A ˆ d [ k − 1] + L 1 ( y [ k ] − ˆ y [ k ]) ˆ y [ k ] = C ˆ x m [ k ] + c l [ k ] + C ˆ d [ k − 1] e [ k ] = y [ k ] − ˆ y [ k ] ˆ d [ k ] = ˆ d [ k − 1] + L 2 C > e [ k ] (11) where ˆ x m [ k ] ∈ R 2 is a state estimate of corrected state vector x m [ k ] at time k ; ˆ y [ k ] ∈ R 2 N is the estimate vector of observ ation y [ k ] ; ˆ d [ k ] is an estimate of the spoofing attack d [ k ] . W e note here that ˆ d [ − 1] , ˆ x m [0] , and ˆ y [0] should be initialized before iteration starts at k = 0 with arbitrary initial conditions. Matrices L 1 ∈ R 2 × 2 N and L 2 ∈ R 2 × 2 are optimization v ariables where L 1 is akin, in principle, to a Luenberger gain that is designed here to ensure robustness of the state estimation to spoofing attacks. B. Design of Robust Gains L 1 , 2 The design of the robust estimator gains L 1 , 2 is based on linear matrix inequalities (LMIs). Simply put, the objective of the designed observ er is to guarantee asymptotically stable estimation error dynamics. That is, matrices L 1 and L 2 are designed to guarantee that lim k →∞ e [ k ] = 0 under non- zero spoofing attack d [ k ] and bounded estimation error un- der spoofing attacks. The R O B U S T E S T I M ATO R variables are designed via solving this low-dimensional feasibility problem with one linear matrix inequality (LMI), gi ven as follows: E S T I M A T O R D E S I G N find G ∈ R 2 × 2 N , P ∈ R 2 × 2 , Q ∈ R 2 N × 2 N , M ∈ R 2 × 2 s . t .     P ? ? ? 0 Q ? ? GC − P A GC − P A P ? M C > C M C > C − Q 0 Q      0 (12a) { P , Q , M } = { P > , Q > , M > }  0 , (12b) where the symbol ? is used to represent symmetric components in symmetric block matrices. After solving ( 12 ) for positive definite matrix variables P , Q , and M , and real matrix variable G , the observ er gains are computed as follows L 1 = P − 1 G , L 2 = M Q − 1 . (13) As mentioned earlier , the state estimator design is derived from [ 19 ]; the reader is referred to that paper for the deriv ation of the above LMIs. Note that no tuning is required to solve E S T I M A T O R D E S I G N , and this LMI can be solved analytically via ev aluating the Karush-Kuhn-T ucker conditions for feasi- bility [ 24 ]. Furthermore, any con ve x optimization toolbox or LMI solv er can be used to solve ( 12 ). These include Matlab’ s LMI solver , CVX [ 25 ], and Y almip [ 26 ]. W e note the follo wing. First, the necessary conditions for existence of L 1 , 2 are standard. These conditions are (a) the detectability of ( A , C ) ; (b) the classical rank matching condi- tion stating that rank( C B d ) = rank( B d ) where B d = I 2 is the matrix coef ficient of d [ k ] in ( 8 ); and (c) bounded variations of the unknown spoofing signal d [ k ] . Second, this robust estimator not only estimates x [ k ] , but also the spoofing attack d [ k ] . This is instrumental in mitigating and correcting the attack. The next section showcases the design of the gain matrices. C. State Corr ection Algorithm under Spoofing Upon solving the LMIs and running the R O B U S T E S T I - M A T O R from any arbitrary initial conditions, the estimator is guaranteed theoretically to produce bounded estimation error e [ k ] , thereby estimating the then-spoofed bias and drift, and reconstructing the spoofing attack ˆ d [ k ] . W ith that in mind, this does not indicate that the bias and drift are correctly 4 IEEE CONTR OL SYSTEMS LETTERS, JUNE 2019, IN PRESS Algorithm 1: Robust Bias and Drift Estimation 1 input: Number of satellites N , matrices A and C 2 initialize: ˆ d [ − 1] , ˆ x m [0] , ˆ y [0] 3 Of fline Computations 4 Compute G, P , Q, M giv en A and C by solving ( 12 ) 5 Obtain L 1 and L 2 from ( 13 ) 6 Online Computations 7 while k ≥ 0 do 8 Obtain c l [ k ] from satellite measurements 9 Run ( 11 ) and obtain ˆ x m [ k + 1] and ˆ d [ k ] 10 Perform spoofing attack correction ( 15 ) 11 Perform state and output correction ( 16 ) 12 k ← k + 1 13 output: Attack-free ˆ x c [ k ] , TSA estimates ˆ d c [ k ] estimated seeing that spoofing attack had already changed the state through the state propagation and difference equation. T o that end, this section develops a state correction algorithm to recov er the authentic bias and drift values in real-time. T o that end, the dynamics ( 9 ) can be written as x m [ k + 1] = Ax m [ k ] + Ad [ k − 1] = A ( Ax m [ k − 1] + Ad [ k − 2]) + Ad [ k − 1] . (14) This relationship rev eals that current spoofing during one time instant comprises the cumulativ e attacks from the pre- vious time steps. Consequently , the disturbance estimate ˆ d is not sufficient enough to correct the attacked states. Rather , a new disturbance estimate vector d c is formulated to account for the accumulated disturbances. Considering the estimate of spoofing attacks ˆ d [ k ] =  ˆ d 1 [ k ] ˆ d 2 [ k ]  > computed by ( 11 ), we propose estimating the new disturbance estimate vector d c via ˆ d c [ k ] =       k X l =1 ˆ d 1 [ l ] + k − 1 X l =1 ( k − l ) ˆ d 2 [ l ] k − 1 X l =1 ˆ d 2 [ l ]       (15) This equation acknowledges the fact that estimated state x m is still contaminated by the attack from the past time step. Therefore, the corrected state x c and authentic observ ation state y c could be retrie ved by subtracting ˆ d c from x m and ˆ y as follows: ˆ x c [ k ] = ˆ x m [ k ] − ˆ d c [ k ] , ˆ y c [ k ] = ˆ y [ k ] − C ˆ d c [ k ] (16) Algorithm 1 showcases the overall problem design, robust state estimation, and the reconstruction of the corrected bias and drift of the GPS recei ver . The algorithm takes as inputs: the fixed number of satellites N , A and C , and satellites data which is encoded through c l [ k ] . The algorithm is di vided into two stages—an offline stage and an online one. In the offline stage, the R O B U S T E S T I M A T O R gains L 1 and L 2 are computed via solving ( 11 ) and ev aluating ( 13 ). The online stage includes running the R O B U S T E S T I M A T O R ( 11 ) and the cumulati ve corrections for the states and spoofing attacks. The algorithm returns the attack-free, yet still slightly noisy T ABLE I V E CT O R N O M EN C L A T U R E . Notation Description x GT ground truth state vector x EKF estimated state v ector from EKF x Luen estimated state v ector from Luenber ger Observer ˆ x m modified state v ector estimates ˆ x c corrected state v ector via ( 16 ) d attack vector vector applied to x ˆ d c corrected spoofing attack vector via ( 15 ) estimates of the bias and drift ˆ x c [ k ] and an estimate of the actual spoofing attack ˆ d c [ k ] . It is note worthy to mention the follo wing. First, the offline component of the algorithm—albeit offline—can be solv ed analytically seeing that the problem dimension is very small, when considering that only fe w satellite measurements are needed. Second, the algorithm and the LMI feasibility problem both only require a fixed number of satellites, rather than a fixed satellite combination. This is important considering that different satellites are visible each time. In short, the proposed algorithm in this paper only assumes a minimum fixed number of satellites N , where these N satellites can be changing in real time without impacting the algorithm or the design of the robust estimator . Third, Algorithm 1 works for any reasonable initial condi- tions, that is, the estimation should con verge regardless of the initial conditions choice. Fourth, this method is truly tuning- free: no prior kno wledge of the statistical distribution of noise, or prior knowledge or tuning of any parameters is needed. The algorithm is also low-memory , as the only computation needed to be performed online is running the R O B U S T E S T I M ATO R and the correction models—both require a small number of matrix-vector multiplications. This implies that the proposed algorithm can be implemented in low-memory de vices without the need for any intensiv e computational effort or internet connectivity . Finally , we note that the proposed algorithm has no stopping criterion seeing that it runs in real time. I V . C A S E S T U D I E S : A R E A L I S T I C T E S T B E D This section discusses the detection and mitigation of TSAs via various approaches. First, the experimental procedure is discussed. Then, we compare the performance of the extended Kalman filter (EKF)—which has long been used in the literature [ 20 ] as a gr ound truth for estimating the bias and drift—and the classical Luenberger observer under spoofing attacks. Then, the performance of the proposed robust estimator under TSAs is showcased, followed by thorough comparison of the performance of the approaches. The following link includes all codes and data used to generate the results, including the acquired GPS data: github.com/junhwanlee95/Robust-Estimator . T able I summarizes the important v ector nomenclature. A. Setup: Model Simulation & Obtaining Raw GPS Data A Google Ne xus 9 tablet, which has an embedded GPS chipset, is used to collect real GPS signals, which are recorded on Nov ember 4, 2018 at the Univ ersity of T exas at San Antonio main campus. The data are available for the reader through the aforementioned Github link. While the receiv er LEE ET AL. ; TUNING-FREE, LOW MEMOR Y R OBUST ESTIMA TOR TO MITIGA TE GPS SPOOFING A TT ACKS 5 acquires the signal, the de vice remained still to simulate the stationary scenario. Raw GPS data is post-processed to obtain pseudorange and pseudorange rate data by GNSS Logger , the Android application released by the Google Android location team [ 27 ]. Then, T ype I and II attacks are injected into the pseudorange and pseudorange rate data to simulate spoofing as discussed in Section II and sho wn in [ 11 , Section III]. The initial conditions for the rob ust estimator are chosen to be different than the actual, ground truth conditions. B. EKF and Luenberger P erformances Under Attacks Here, we are interested in testing whether the EKF and the classical Luenberger observer [ 28 ] can withstand TSAs of T ype I and II. After T ype I and II attacks are applied from t = 30 s to t = 400 s , the performance of the EKF and Luenberger observer—which does not assume any statistical distribution about the noise—are shown in Fig. 1 and 2 . While the ground truth clock bias and drift, x GT , are ac- quired through EKF by processing the authentic pseudorange data, the x EKF bias and drift are generated by applying the EKF to the spoofed pseudoranges. Another comparison to x GT is offered by the Luenberger observer estimate x Luen produced after designing the observer gain L 1 such that the closed loop system eigen values are at 0 . 5 and 0 . 7 . The performances of EKF and Luenberger observer are shown respectiv ely in Fig. 1 and Fig. 2 . It is evident that both approaches fail to estimate the correct states in the presence of the attack. C. Rob ust Estimator P erformance under T ype I/II TSAs In this section, we test the proposed robust estimator and run Algorithm 1 which contains an offline stage and an online one. First, we set N = 4 , i.e., we choose to sample data from only four satellites. W e solve the LMIs ( 12 ) for the estimator gains. Using these estimator gains, the online portion of Algorithm 1 is run. Pertaining to T ype I attack, Fig. 3 showcases the performance of R O BU S T E S T I M ATO R . Due to the attack at t = 30 s, the ˆ x c clock estimates are initially not correct, but the clock bias and drift approach the respectiv e ground truth values within approximately 3 and 11 seconds (cf. Fig. 3 (a) and 3 (b) respectiv ely). Under the same condition and procedure, Algorithm 1 is applied on T ype II attack to detect and correct the spoofed states. The results of obtaining the corrected state ˆ x c are sho wn in Fig. 4 . In order to accurately depict the performance of R O - B U S T E S T I M A T O R , the relati ve estimation error is calculated as | ˆ x c − x GT | x GT for each of the two states in x ov er time. The resulting graphs are sho wn in Fig. 5 . Comparison between ˆ x c and the x GT rev eals that the maximum error between the two biases is 952 . 09 m or 3 . 17 µ s . It is thus demonstrated that the Robust Estimator of Algorithm 1 successfully detects and corrects the T ype II attack. In the interest of gauging the performance of each approach, the root mean square error ( RMSE ) of the estimated clock bias is calculated under both attack types. Let K denote the total length of observation time ( K = 400 in this experiment). The RMSE is defined as RMSE = q 1 K P K − 1 k =0 ( c ˆ b u [ k ] − c ˇ b u [ k ]) 2 0 100 200 300 400 Time(s) (a) -4 -2 0 10 4 0 100 200 300 400 Time(s) (b) -10000 -5000 0 Fig. 1. Comparison of ground truth, EKF and Luenberger observer estimates under T ype I attack: (a) clock bias; (b) clock drift. 0 50 100 150 200 250 300 350 400 Time(s) (a) -2 -1 0 1 2 10 5 0 100 200 300 400 Time(s) (b) -1000 -500 0 500 1000 Fig. 2. Comparison of ground truth, EKF and Luenberger observer estimates under T ype II attack: (a) clock bias; (b) clock drift. where c ˇ b u [ k ] is the ground truth clock bias under normal conditions, and c ˆ b u [ k ] equals to the estimated clock bias value from each approach. Under T ype II attack, the RMSE for EKF and the Luenberger observer are RMSE EKF = 74344 m and RMSE Luenberger = 74433 m respectiv ely , while that of the robust estimator is RMSE RobustEstimator = 354 . 9 m . As for T ype I attack, the RMSEs are as follo ws: RMSE EKF = 8477 . 3 m , RMSE Luenberger = 8104 . 9 m , and RMSE RobustEstimator = 1029 m . This illustrates the perfor- mance of this tuning-free, lo w-memory rob ust estimator in detecting spoofing attacks, while correctly reconstructing the bias and drift states of the GPS recei ver . V . P A P E R S U M M A RY A N D F U T U R E W O R K In this paper , the design and realistic application of a low- memory , real-time R O BU S T E S T I M AT O R is studied. Utilizing the GPS receiv er on a Google Nexus 9, real GPS data are collected and post-processed by injecting time-synchronization attacks to spoof the clock bias and drift of the device. T wo types of attacks are introduced, and tested by the designed estimator . The estimator successfully detects and estimates the spoofing attacks on each state, and mitigates the spoofing on both types of attack by furnishing the corrected clock states to 6 IEEE CONTR OL SYSTEMS LETTERS, JUNE 2019, IN PRESS 0 50 100 150 200 250 300 350 400 (a) -4 -2 0 Bias (m) 10 4 350 400 -4 -3.5 10 4 0 50 100 150 200 250 300 350 400 Time (s) (b) -100 -90 -80 Drift (m/s) Fig. 3. Comparison of ground truth and corrected state through the correc- tion ( 16 ) under T ype I attack: (a) clock bias; (b) clock drift. 0 50 100 150 200 250 300 350 400 (a) -4 -2 0 Bias (m) 10 4 350 400 -4 -3.5 10 4 0 50 100 150 200 250 300 350 400 Time (s) (b) -100 -90 -80 Drift (m/s) Fig. 4. Comparison of ground truth and corrected state through the correc- tion ( 16 ) under T ype II attack: (a) clock bias; (b) clock drift. the user . Future work will focus on developing robust estima- tors under spoofing attacks for non-stationary GPS recei vers, which inv olve nonlinearities in the GPS measurement model. R E F E R E N C E S [1] U.S. Gov ernment Accountability Office, “GPS disruptions: Efforts to assess risks to critical infrastructure and coordinate agency actions should be enhanced, ” GA O-14-15, Nov . 2014. [Online]. A v ailable: https://www .gao.gov/products/GA O- 14- 15 [2] X. Li and Q. Xu, “ A reliable fusion positioning strategy for land vehicles in GPS-denied en vironments based on lo w-cost sensors, ” IEEE T rans. Ind. Electr on. , vol. 64, no. 4, pp. 3205–3215, Apr . 2017. [3] D. Schmidt, K. Radke, S. Camtepe, E. Foo, and M. Ren, “ A survey and analysis of the GNSS spoofing threat and countermeasures, ” ACM Computing Surve ys , v ol. 48, no. 4, pp. 1–31, May 2016. [4] T . Nighswander , B. Ledvina, J. Diamond, R. Brumley , and D. Brumley , “GPS software attacks, ” in Pr oc. of the A CM Conf. on Comput. and Commun. Security , Raleigh, NC, Oct. 2012, pp. 450–461. [5] M. L. Psiaki and T . E. Humphreys, “GNSS spoofing and detection, ” Pr oceedings of the IEEE , vol. 104, no. 6, pp. 1258–1270, June 2016. [6] L. Heng, J. J. Makela, A. D. Domnguez-Garca, R. B. Bobba, W . H. Sanders, and G. X. Gao, “Reliable gps-based timing for power systems: A multi-layered multi-receiver architecture, ” in Proc. P ower and Ener gy Confer ence at Illinois (PECI) , Champaign, IL, Feb .-Mar. 2014, pp. 1–7. [7] T . E. Humphreys, B. M. Ledvina, M. Psiaki, B. W . OHanlon, and J. P . M. Kintner , “ Assessing the spoofing threat: Development of a portable GPS civilian spoofer, ” in Pr oc. 21st Int. T echnical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS) , Sav anna, GA, Sept. 2008, pp. 2314–2325. [8] B. Moussa, M. Debbabi, and C. Assi, “Security assessment of time synchronization mechanisms for the smart grid, ” IEEE Commun. Surve ys T ut. , vol. 18, no. 3, pp. 1952–1973, thirdquarter 2016. [9] F . Zhu, A. Y oussef, and W . Hamouda, “Detection techniques for data- lev el spoofing in GPS-based phasor measurement units, ” in Pr oc. Int. Conf. on Selected T opics in Mobile Wir eless Networking (MoWNeT) , Cairo, Egypt, Apr . 2016, pp. 1–8. 0 50 100 150 200 250 300 350 400 Time(s) 0 1 2 3 4 Relative Error 10 -3 Fig. 5. Relative estimation error for clock bias and drift defined as: RE 1 = | ˆ x c (1) − x GT (1) | x GT (1) , RE 2 = | ˆ x c (2) − x GT (2) | x GT (2) . [10] S. Han, D. Luo, W . Meng, and C. Li, “ A novel anti-spoofing method based on particle filter for GNSS, ” in Pr oc. IEEE Int. Conf. on Com- munications (ICC) , Sydne y , Australia, June 2014, pp. 5413–5418. [11] A. Khalajmehrabadi, N. Gatsis, D. Akopian, and A. F . T aha, “Real- time rejection and mitigation of time synchronization attacks on the global positioning system, ” IEEE Tr ansactions on Industrial Electronics , vol. 65, no. 8, pp. 6425–6435, Aug. 2018. [12] D.-Y . Y u, A. Ranganathan, T . Locher , S. ˇ Capkun, and D. Basin, “Short paper: Detection of GPS spoofing attacks in po wer grids, ” in Proc. of the ACM Conf. on Security and Privacy in W ireless Mobile Networks , Oxford, United Kingdom, July 2014, pp. 99–104. [13] K. Jansen, N. O. Tippenhauer , and C. P ¨ opper , “Multi-receiv er GPS spoofing detection: Error models and realization, ” in Proc. 32nd Annual Confer ence on Computer Security Applications , Los Angeles, CA, Dec. 2016, pp. 237–250. [14] S. Bhattacharyya, “Observer design for linear systems with unknown inputs, ” IEEE T rans. Autom. Control , vol. 23, no. 3, pp. 483–484, June 1978. [15] W . Kang, A. J. Krener, M. Xiao, and L. Xu, “ A survey of observers for nonlinear dynamical systems, ” in Data Assimilation for Atmospheric, Oceanic and Hydr ologic Applications (V ol. II) , S. Park and L. Xu, Eds. Berlin, Germany: Springer, 2013, pp. 1–25. [16] C. Masreliez and R. Martin, “Rob ust Bayesian estimation for the linear model and rob ustifying the Kalman filter, ” IEEE Tr ans. Autom. Contr ol , vol. 22, no. 3, pp. 361–371, June 1977. [17] A. H. Jazwinski, Stoc hastic Pr ocesses and Filtering Theory . Mineola, NY : Dover Publications, 2007. [18] G. A. Einicke and L. B. White, “Robust extended kalman filtering, ” IEEE T rans. Signal Processing , vol. 47, no. 9, pp. 2596–2599, Sept. 1999. [19] B. A. Charandabi and H. J. Marquez, “Observer design for discrete-time linear systems with unknown disturbances, ” in Pr oc. IEEE 51st IEEE Conf. Decision and Contr ol (CDC) , Maui, HI, Dec. 2012, pp. 2563– 2568. [20] P . Axelrad and R. G. Brown, “GPS na vigation algorithms, ” in Global P ositioning System: Theory and Applications , B. W . Parkinson, J. J. Spilker , P . Axelrad, and P . Enge, Eds. American Institute of Aeronautics and Astronautics, 1996, vol. I, ch. 9. [21] D. P . Shepard, T . E. Humphreys, and A. A. Fansler , “Ev aluation of the vulnerability of phasor measurement units to GPS spoofing attacks, ” Int. J. Crit. Infrastruct. Pr otect. , vol. 5, no. 3, pp. 146–153, 2012. [22] T . Humphreys, J. Bhatti, D. Shepard, and K. W esson, “The Texas spoofing test battery: T oward a standard for ev aluating GPS signal authentication techniques, ” in Proc. 25th Int. T echnical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS) , Nashville, TN, Sept. 2012, pp. 3569–3583. [23] IEEE Standar d for Synchr ophasor Measur ements for P ower Systems , IEEE Std. C37.118.1-2011 (Re vision of IEEE Std. C37.118-2005), 2011. [24] S. Boyd and L. V andenber ghe, Conve x Optimization . Cambridge Univ ersity Press, 2004. [25] M. Grant and S. Boyd, “CVX: Matlab software for disciplined con vex programming, version 2.1, ” http://cvxr .com/cvx , 2014. [26] J. L ¨ ofberg, “Y almip: A toolbox for modeling and optimization in matlab, ” in Proceedings of the CA CSD Confer ence , vol. 3. T aipei, T aiwan, 2004. [27] “ Android GNSS, ” https://developer .android.com/guide/topics/sensors/ gnss.html , accessed: 2017-02-20. [28] D. Luenberger , “Observers for multivariable systems, ” IEEE Tr ans. Autom. Control , vol. 11, no. 2, pp. 190–197, Apr . 1966.