GPS Spoofing Mitigation and Timing Risk Analysis in Networked PMUs via Stochastic Reachability
GGPS Spoofing Mitigation and Timing RiskAnalysis in Networked PMUs via StochasticReachability
Sriramya Bhamidipati,
University of Illinois at Urbana-Champaign
Grace Xingxin Gao,
Stanford University
BIOGRAPHIESSriramya Bhamidipati is a Ph.D. student in the Department of Aerospace Engineering at the University of Illinois at Urbana-Champaign, where she also received her master’s degree in 2017. She obtained her B.Tech. in Aerospace from the IndianInstitute of Technology, Bombay in 2015. Her research interests include GPS, power and control systems, artificial intelli-gence, computer vision, and unmanned aerial vehicles.
Grace Xingxin Gao is an assistant professor in the Department of Aeronautics and Astronautics at Stanford University. Beforejoining Stanford University, she was an assistant professor at University of Illinois at Urbana-Champaign. She obtained herPh.D. degree at Stanford University. Her research is on robust and secure positioning, navigation and timing with applicationsto manned and unmanned aerial vehicles, robotics, and power systems.
ABSTRACT
To address PMU vulnerability against spoofing, we propose a set-valued state estimation technique known as StochasticReachability-based Distributed Kalman Filter (SR-DKF) that computes secure GPS timing across a network of receivers. Uti-lizing stochastic reachability, we estimate not only GPS time but also its stochastic reachable set, which is parameterized viaprobabilistic zonotope (p-Zonotope). While requiring known measurement error bounds in only non-spoofed conditions, wedesign a two-tier approach: We first perform measurement-level spoofing mitigation via deviation of measurement innovationfrom its expected p-Zonotope, and second perform state-level timing risk analysis via intersection probability of estimated p-Zonotope with an unsafe set that violates IEEE C37.118.1a-2014 standards. We validate the proposed SR-DKF by subjectinga simulated receiver network to coordinated signal-level spoofing. We demonstrate an improved GPS timing accuracy and suc-cessful spoofing mitigation via our SR-DKF. We validate the robustness of the estimated timing risk as the number of receiversare varied.
Modern power systems analyze real-time grid stability using a widely dispersed network of advanced devices, known as PhasorMeasurement Units (PMUs) [1]. As per the American Recovery and Reinvestment Act of , almost
PMUs have beendeployed across the US and Canada as of March [2, 3]. Figure 1 shows a map of the dense network of these spatiallydistributed PMUs. The North American Synchrophasor Initiative network (NASPInet), which is a collaborative effort by theU.S. Department of Energy and the North American Electric Reliability Corporation, provides a secure and decentralizedinfrastructure for communication among the PMUs [4]. Each PMU is equipped with a GPS receiver to precisely time-stamp thevoltage and current phasor values at critical power substations. A network of geographically distributed GPS receivers facilitaterobust time authentication by cross-checking each other and thereby eliminate the need for an external reference receiver [5].The networked approach in power systems overlays the existing grid infrastructure and therefore can be implemented withoutadditional hardware.PMUs are susceptible to GPS spoofing because of the low signal power and unencrypted structure of the civilian GPS/L1signals [6]. As validated in prior literature [7], this vulnerability can be exploited by attackers to mislead the GPS timing at one a r X i v : . [ c s . M A ] J a n r more critical substations and cause cascading system failure. The various spoofing attacks that can potentially disrupt GPStiming vary from a simple meaconing to more sophisticated signal-level spoofing and coordinated attacks. A spoofer executingmeaconing (or time jump) broadcasts earlier recorded authentic GPS signals to induce a constant lag in the victim receiver’scomputed time as compared to the true time. In signal-level spoofing (or time walk), the spoofer broadcasts the simulated GPSsignals to gradually steer the computed time of victim receiver away from its true time. In addition, coordinated attack occurswhen multiple spoofers manipulate the GPS signals received at one or more victim receivers. We refer readers to our prior workfor a detailed overview of these attacks [8].Secure GPS timing is of paramount importance for reliable PMU operations. The IEEE C37.118.1a-2014 standards forSynchrophasors [9] define % Total Vector Error (TVE) as a benchmark for grid stability. For a PMU operating at Hz,which is a commonly used frequency for phasor data logging, > % TVE is caused by a timing error of > ± . µ s. Thekey measures of secure GPS timing are spoofing attack mitigation and timing risk analysis . In this context, we define risk asthe probability of estimation error to exceed a safety Alert Limit (AL). In this work, we set AL = ± . µ s to ensure that theestimation error in GPS time complies with the IEEE C37.118.1a-2014 standards.Figure 1: A network of > PMUs deployed across the US [10]The state-of-the-art prior work [11–13] on networked anti-spoofing falls under the category of point-valued state estimationapproaches. A point-valued state estimate denotes the expected value of the state vector given the measurements. The existingliterature on the point-valued anti-spoofing methods, which is summarized in [8, 14, 15], includes signal quality monitoring,angle-of-arrival and time-difference-of-arrival-based analyses of the incoming satellite signals, validation of the encryptioncodes, and attack magnitude estimation. Unlike the point-valued techniques, a set-valued approach utilizes the pre-computed(or known) error bounds on system dynamics and received measurements to propagate the set of state estimates [16]. Anillustration of the difference between point-valued and set-valued propagation is seen in Fig. 2(a). (a) Set-valued approach (b) Stochastic reachability
Figure 2: Set-valued approach and stochastic reachability; (a) shows an illustration of the difference between point-valued andset-valued approach. Unlike the point-valued techniques where only the point-valued state is propagated across time, the set-valued approach propagates the set of state estimates, thereby providing confidence bounds on the estimation errors; (b) showsan intuitive understanding of stochastic reachability that computes a set of stochastic reachable states given an initial stochasticset and known error bounds.Much prior literature is available on set-valued theory [16–18]. Our current work is inspired by a set-valued formulation2esigned for stochastic processes known as stochastic reachability [19]. Stochastic reachability is widely used in robotics forpath planning and collision avoidance, however, it has not been applied in the field of power systems. As shown in Fig. 2(b),stochastic reachability computes a set of stochastic reachable states given an initial stochastic set of states, and the knownstochastic sets of system disturbances and measurement errors. In this context, each state in the stochastic reachable set isassociated with a probabilistic measure that indicates the likelihood of its occurrence. Our choice of stochastic reachability isjustified because it accounts for the stochastic errors in receiver clock dynamics, data transfer latencies, and GPS measurements.In this work, we propose a set-valued state estimation technique for a network of GPS receivers, termed as StochasticReachability-based Distributed Kalman Filter (SR-DKF). The stochastic reachability formulation requires known error bounds,however, the attack errors induced via spoofing are unbounded and time-varying in nature. To account for this, we consider theerror bounds in GPS measurements to be known only in non-spoofed (or authentic) conditions and pre-compute these boundsfrom offline empirical analysis of historical data [20]. This current work is based on our recent ION GNSS+ conferencepaper [21].In a conventional DKF approach [22], each receiver evaluates the received measurements against a point-valued attackthreshold and thereafter fuses the measurements obtained from different receivers to compute its point-valued state estimate. Incontrast, our SR-DKF performs secure state estimation of GPS time and spoofing attack mitigation in the set-valued domain.Therefore, we estimate not only the point-valued mean but also the stochastic reachable set of GPS time that is later utilizedto compute the associated risk of violating IEEE C37.118.1a-2014 standards. Unlike the point-valued state estimation tech-niques [8,12], our proposed SR-DKF algorithm does not require the conservative assumption that all satellite measurements at avictim receiver are attacked by the spoofer. Our set-valued analysis of spoofing attack combined with the redundancy in numberof available measurements from a network of GPS receivers, eliminates the need to estimate the attack magnitude [23]. To ourknowledge, no prior literature exists on performing the secure state estimation using stochastic reachability nor on quantifyingthe timing risk associated with estimated GPS time.
The main contributions of this paper are listed as follows:1. Utilizing an efficient set representation known as probabilistic zonotope (p-Zonotope) [24] to parameterize the stochasticreachable set, we develop a set-valued DKF framework to compute secure GPS time for a network of receivers.2. Considering offline-estimated measurement error bounds in non-spoofed (or authentic) conditions, we design a two-tierapproach that performs(a)
Measurement-level spoofing attack mitigation by adaptively evaluating the deviation of point-valued GPS mea-surements from the set-valued domain of measurement innovation. Therefore, the proposed SR-DKF algorithm isresilient to a variety of GPS attacks, such as meaconing, signal-level spoofing, coordinated attack, etc.(b)
State-level timing risk analysis by estimating the intersection probability of the estimated set of stochastic reachablestates (represented by p-Zonotope) with an unsafe set, i.e., a set violating the IEEE C37.118.1a-2014 standards [9].3. For a simulated setup, we demonstrate the successful mitigation of both simple meaconing and sophisticated signal-levelspoofing attacks. During a simulated case of coordinated attack that affects multiple receivers, we show an improvementin the state estimation accuracy via the proposed SR-DKF algorithm as compared to conventional point-valued DKF andsingle-receiver-based adaptive KF. We also demonstrate the robustness of the proposed SR-DKF algorithm by analyzingthe variation in timing risk with respect to the number of receivers and magnitude of spoofing attacks.The rest of the paper is organized as follows: Section 2 outlines the system model, the preliminaries of stochastic reacha-bility and the conventional point-valued DKF; Section 3 describes the proposed SR-DKF algorithm; Section 4 experimentallyvalidates the improved state estimation accuracy of the GPS timing computed via SR-DKF as compared to the point-valuedDKF, and also demonstrates the trend of timing risk associated with the estimated GPS time; and Section 5 concludes the paper.
The three main design aspects of the proposed SR-DKF algorithm, namely i) a spatially distributed network of static GPSreceivers, ii) a distributed processing framework, and iii) the set representation via p-Zonotopes, are described as follows:3 ) A spatially distributed network of static GPS receivers:
We utilize the existing grid infrastructure to consider a handful (around four to eight) of the spatially distributed PMUs (orGPS receivers). This framework provides two advantages: First, by considering a spatially distributed network, we increase themeasurement redundancy while minimizing the probability of multiple receivers being simultaneously spoofed; Second, giventhat the power substations are static, we can utilize the known receiver positions to aid the proposed SR-DKF algorithm. ii) A distributed processing framework:
We opt for distributed processing to minimize the reaction time to attacks and also to localize the extent of failure. This isjustified because global GPS reference time for synchronizing the PMUs is no longer reliable during spoofing. Therefore, theconventional communication protocol involving send or receive requests cannot be scheduled across the network. Each receiverin our distributed network asynchronously broadcasts the latest data to its neighbors, thereby reducing communication latency. (a) Example: 1D p-Zonotope (b) Example: 2D p-Zonotope
Figure 3: Illustration of p-Zonotopes; (a) shows the example of a 1D p-Zonotope and is indicated by the Parula colormap [25].This 1D p-Zonotope encompasses multiple Gaussian distributions, a subset of which are shown in gray; (b) shows the exampleof a 2D p-Zonotope and is indicated by the Parula colormap. This 2D p-Zonotope is used to represent the state vector of interestto us in this paper and whose details are given later in Section 3. iii) The set representation via p-Zonotopes:
We represent stochastic reachable set via an enclosing probabilistic hull [24] known as p-Zonotope, denoted by L . Amongall available set representations [26] such as zonotope, ellipsoid, polytope, etc., our choice of p-Zonotope is justified becauseit can efficiently enclose the stochastic biases in GPS measurements. A p-Zonotope can enclose a variety of distributionssuch as mixture models, distributions with uncertain/time-dependent mean, non-Gaussian/unknown distribution, etc. A one-Dimensional (1D) example of a p-Zonotope that encloses multiple Gaussian distributions with varying mean and covariance isillustrated in Fig. 3(a). The probabilistic nature of p-Zonotopes blends well with the Gaussian-centric framework of DKF.Unlike the conventional techniques that propagate point-valued state estimates, the proposed SR-DKF algorithm utilizesthe pre-computed error bounds on the receiver clock dynamics and measurements to recursively propagate the p-Zonotope ofthe state estimate. An example of a two-Dimensional (2D) p-Zonotope is shown in Fig. 3(b). In this work, the 2D p-Zonotopeis used to represent the state estimate that is introduced later in Section 3. A p-Zonotope L , as seen in Eq. (1a) and Fig. 3(a), is characterized by three parameters [24], namely, a) the mean of its center,denoted by c ∈ R n , b) the uncertainty in the center, represented by the generator matrix G ∈ R n × e , and c) the over-boundingcovariance of the Gaussian tails, denoted by Σ ∈ R n × n . L = ( c , G, Σ) (1a)4ntuitively, a p-Zonotope with a certain center mean, i.e., zero center uncertainty G = , represents a Gaussian distribution thatoverbounds multiple distributions with the same center mean c . This is mathematically represented by Z in Eq. (1b). Z = (cid:40) c + l (cid:88) i =1 N i (0 , g i (cid:41) , (1b)where N i (0 , ∀ i ∈ { , · · · , l } denotes independent normally distributed random variables and g i ∈ R n denotes the associatedgenerator vectors of Z , such that G = [ g , · · · , g l ] and Σ = GG (cid:62) .Similarly, a p-Zonotope with zero covariance Σ = simplifies to a zonotope [26] that encompasses all the possible valuesof the center, such that the generator matrix of Z is G = [ g , · · · , g e ] . This case is mathematically represented by Z in Eq. (1c). Z = (cid:68) c , G (cid:69) = (cid:40) c + e (cid:88) i =1 β i g i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ≤ β i ≤ (cid:41) (1c)The p-Zonotope is a combination of both sets, i.e., L = Z (cid:1) Z , where (cid:1) is a set-addition operator, such that L = sup (cid:8) f Z (cid:12)(cid:12)(cid:12) E [ Z ] ∈ Z (cid:9) , where f Z is the Probability Density Function (PDF) of the distributions represented by Z and E [ · ] denotes the expectation operator. From Eqs. (1a)-(1c), note that L depends on both G and G . Additionally, note that unlikea PDF, the area enclosed by a p-Zonotope does not equal to one. More details regarding the characteristics of zonotopes andp-Zonotopes can be found in the prior literature [23, 24].For example, the 2D p-Zonotope shown in Fig. 3(b) is formulated by considering the bounds on the mean and covariance ofthe first state to be [ − , and [0 , , respectively, and of the second state to be [ − , and [0 , , respectively. Mathematically,we represent this 2D p-Zonotope L = ( c , G, Σ) by c = (cid:20) (cid:21) , G = (cid:20) (cid:21) and Σ = 3 × (cid:20) (cid:21) = (cid:20) (cid:21) . We consider amultiplication factor in the over-bounding covariance Σ based on the empirical rule according to which . of the valueslie within three standard deviations of the mean [27].To perform set operations on the stochastic reachable sets, we require three important properties [24]:• Minkowski sum of addition:
The addition of two p-Zonotopes L = ( c , G , Σ ) and L = ( c , G , Σ ) is given by L (cid:77) L = ( c + c , [ G , G ] , Σ + Σ ) , (2a) p (cid:77) i =1 L i = ( p (cid:88) i =1 c i , [ G , · · · , G p ] , p (cid:88) i =1 Σ i ) , (2b)where (cid:76) is known as the Minkowski sum operator and [ · , · ] denotes a horizontal concatenation operator.• Linear map:
Given a matrix A ∈ R n × n and a p-Zonotope L = ( c , G, Σ) , A L = ( A c , AG, A Σ A (cid:62) ) (2c)• Translation:
Given a vector µ ∈ R n and a p-Zonotope L = ( c , G, Σ) , µ + L = ( µ + c , G, Σ) (2d) Based on the key features of the SR-DKF algorithm, we mathematically represent the inter-connected network of power sub-stations, each mounted with one GPS receiver, as a undirected graph G = { M, E } , where M indicates the set of nodes and E denotes the set of connections between the nodes. Each node represents a GPS receiver in the network and each edge indicatesthe presence of a secure bidirectional communication link between the power substations. The total number of receivers in thenetwork are (cid:12)(cid:12) M (cid:12)(cid:12) , where (cid:12)(cid:12) · (cid:12)(cid:12) denotes the cardinality of a set. Our SR-DKF algorithm does not require the network of powersubstations to be fully connected, and hence, (cid:12)(cid:12) E (cid:12)(cid:12) ≤ (cid:12)(cid:12) M (cid:12)(cid:12)(cid:0)(cid:12)(cid:12) M (cid:12)(cid:12) − (cid:1) .5he neighborhood set of any i th receiver ∀ i ∈ M is represented by M i and indicates the subset of nodes among the M nodes that can communicate with the i th receiver, including itself. Therefore, the associated cardinality of the neighborhood setis denoted by (cid:12)(cid:12) M i (cid:12)(cid:12) where (cid:12)(cid:12) M i (cid:12)(cid:12) ≥ ∀ i ∈ M . Every receiver in the network has a prior knowledge regarding the pre-computedstatic position of all the receivers in its neighborhood set.We define the global state vector at the k th time iteration, which is represented by x k ∈ R , to be comprised of GPS timeand its drift rate, which are given by T k and ˙ T k , respectively. As explained earlier in Section 1, estimated GPS time T k isused in PMUs for time-tagging the phasor measurements. At the k th iteration, the global state vector is the same across allthe receivers in the network. Additionally, we define the local state vector at the i th receiver by y i,k ∈ R . The local statevector comprises two states, namely, the receiver clock bias and clock drift, which are denoted by cδt k and cδ ˙ t k , respectively.At the i th receiver, the global and local state vectors are related via secure time-varying local clock parameters, i.e., T irx,k and ∆ T irx,k , which denote the measured local clock time and the sampling time, respectively. Therefore, x k = y i,k + θ i,k where θ i,k = [ T irx,k , ∆ T irx,k ] (cid:62) . At the k th iteration, each i th receiver maintains an estimate of both the local and global state vectors.When a receiver is spoofed, the local state vector y i,k is manipulated, thereby misleading the global GPS time T k .The state-space model representing the true system at the i th receiver is given by x k = F x k − + ν i,k , z i,k = H i x k + ω i,k , ∀ i ∈ M, (3)where– z i,k denotes the N i × translated GPS measurement vector, such that z i,k = (cid:104) ˜ ρ i,k , ˜ φ i,k , · · · , ˜ ρ N i i,k , ˜ φ N i i,k (cid:105) , where N i denotes the number of visible GPS satellites from the i th receiver. Also, ˜ ρ ji,k and ˜ φ ji,k denote the GPS pseudorange andDoppler residual associated with the j th satellite ∀ j ∈ { , · · · , N i } , respectively, such that (cid:34) ˜ ρ ji,k ˜ φ ji,k (cid:35) = (cid:34) ρ ji,k φ ji,k (cid:35) + θ i,k − (cid:20) || p i − p j |||| ˙ p i − ˙ p j || (cid:21) . Here, ρ ji,k and φ ji,k denote the measured GPS pseudorange and Doppler value.Also, p i and ˙ p i denote the known position and zero velocity of the i th static receiver, respectively, and p j and ˙ p j denotethe pre-computed position and velocity of the j th satellite, respectively;– F denotes the first-order linear state transition model, such that F = (cid:20) T irx,k (cid:21) , f L is the frequency of the GPS L1signal, and c is the speed of light;– H i denotes the GPS measurement model, such that H i = N i × where a × b denotes the matrix of ones with a rowsand b columns;– ν i,k denotes the associated × process noise vector and ω i,k denotes the N i × measurement noise vector that isassociated with the GPS pseudoranges and Doppler values. Based on [22], the general framework of a point-valued DKF implementation performed independently at the i th receiver, isshown in Eqs. (4)-(6). During initialization, i.e., k < , we define the predicted value of the point-valued mean and covariancethat are given by ˆ x i, − and ˆ P i, − , respectively. The process and measurement noise are modeled via a zero-mean Gaussiandistribution that are given by ν i,k ∼ N ( , Q i,k ) and ω i,k ∼ N ( , R i,k ) , respectively, where Q i,k and R i,k represent the processand measurement noise covariance matrix, respectively. After initialization, the conventional point-valued DKF recursivelyexecutes two steps: one is the measurement update and the other is the time update. ˆ x i, − = E (cid:2) x (cid:3) ˆ P i, − = E (cid:104)(cid:16) x − E (cid:2) x (cid:3)(cid:17)(cid:16) x − E (cid:2) x (cid:3)(cid:17) (cid:62) (cid:105) (4)During the time update, as seen in Eq. (5), the corrected variables ¯ x i,k − and ¯ P i,k − from the previous time, i.e., at the ( k − th iteration, are propagated forward-in-time to predict the estimates for the next k th iteration, which are given by ˆ x i,k and ˆ P i,k , respectively. ˆ x i,k = F ¯ x i,k − ˆ P i,k = F ¯ P i,k − F (cid:62) + Q i,k (5)6t the k th iteration, each i th receiver receives the GPS measurements z j,k and measurement noise covariance matrix R j,k from the other receivers in the neighborhood set, i.e., j ∈ M i , and processes them sequentially during the measurement update.In the measurement update, the expected measurements (obtained via the predicted estimates) and the received measurementsare weighted via an estimated Kalman gain to obtain the corrected estimates of the point-valued mean and covariance of the statevector, i.e., ¯ x i,k and ¯ P i,k , respectively. The optimal Kalman gain denoted by K j,k ∀ j ∈ M i is seen in Eq. (6a). For an adaptiveimplementation of point-valued DKF, additional parameters, namely forgetting factor ψ and pre-measurement residuals (cid:15) i,k ,are considered in Eq. (6b) to adaptively estimate the measurement noise covariance matrix R i,k at each time iteration. ¯ P − i,k = ˆ P − i,k + (cid:88) j ∈ M i (cid:16) H (cid:62) j R − j,k H j (cid:17) , ¯ x i,k = ˆ x i,k + (cid:88) j ∈ M i K j,k (cid:16) z j,k − H j ˆ x i,k (cid:17) , (6a)where (cid:15) j,k = z j,k − H j ˆ x j,k ,R j,k = ψR j,k − + (1 − ψ )( (cid:15) j,k (cid:15) (cid:62) j,k + H j ˆ P j,k H (cid:62) j ) ,K j,k = ¯ P i,k H (cid:62) j R − j,k . (6b) Based on the design aspects discussed in Section 2, we first outline the architecture of the SR-DKF and later describe thestep-by-step algorithmic details.Fig. 4 shows the flowchart of the proposed SR-DKF algorithm that is executed independently at the i th receiver/powersubstation, and comprises four main modules, namely time update of set-valued DKF, spoofing attack mitigation, measurementupdate of set-valued DKF, and timing risk analysis. First, we perform a set-valued time update to compute the predicted p-Zonotope of the global state vector. Next, at each i th receiver, we independently utilize the predicted p-Zonotope and receivedGPS measurements to adaptively estimate its attack status. Thereafter, the GPS measurements and their estimated attack statusare asynchronously broadcast across receivers within the network. In the measurement update, the GPS measurements andtheir attack status from neighbors, i.e., j ∈ M i , are collectively processed to compute the corrected p-Zonotope. The centermean of the corrected p-Zonotope is the estimated point-valued GPS time that is given to PMUs. The center uncertainty andover-bounding covariance are further analyzed for their intersection probability with an unsafe set to compute the associatedtiming risk. In this context, we define unsafe set by a set of global states that violate the IEEE C37.118.1a-2014 standards.Thereafter, the process repeats for the next time iteration.As explained earlier in Section 1 and from prior literature [20], it is justified to consider the error bounds to be knownin the non-spoofed (or authentic) conditions. We can empirically calculate the measurement bounds in authentic conditionsby analyzing the large amounts of historical GPS data from publicly-accessible websites [28], namely Crustal Dynamics DataInformation System (CDDIS), International GNSS Service (IGS), Continuously Operating Reference Station (CORS), etc.Similarly, we can obtain empirical bounds on process noise by evaluating the accuracy of protocols used in power systems forsynchronization across power substations [29] , such as Precise Time Protocol (PTP), Network Time Protocol (NTP), etc.From Eq. (1a), we define the p-Zonotopes of process and measurement noise bounds, which are represented by L i ν and L i ω ,respectively, as follows: L i ν = (cid:0) , G i ν , Σ i ν (cid:1) L i ω = (cid:0) , G i ω , Σ i ω (cid:1) (7)where G i ν and Σ i ν denote the generator and covariance matrices of p-Zonotope that represents process noise bounds. Similarly, G i ω and Σ i ω denote the generator and covariance matrices of p-Zonotope that represents measurement noise bounds. Note thatthe over-bounding covariance matrices Σ i ν and Σ i ω represent the characteristics of p-Zonotopes and are different from Q i,k and R i,k , defined earlier in Section 2.3. 7igure 4: Architecture of the proposed SR-DKF algorithm, which is executed independently at the i th receiver, comprises fourmain modules, namely time update of set-valued DKF, spoofing attack mitigation, measurement update of set-valued DKF, andtiming risk analysis. We formulate the set-valued DKF by applying stochastic reachability to the point-valued DKF described earlier inEqs. (4)-(5). Similar to point-valued DKF, the set-valued DKF also performs time and measurement updates to computethe predicted and corrected stochastic reachable sets of the global state vector x k , respectively. To perform spoofing attackmitigation and timing risk analysis across a network of GPS receivers, we analyze the predicted and corrected stochasticreachable sets of the state estimation error, respectively. At the k th iteration, the estimation error in the point-valued predictedand corrected state estimate is given by ∆ ˆ x i,k and ∆ ¯ x i,k , respectively, where ∆ ˆ x i,k = ˆ x i,k − x k and ∆ ¯ x i,k = ¯ x i,k − x k .We represent the predicted and corrected stochastic reachable sets of state estimation error via predicted and corrected p-Zonotopes, respectively. The predicted and corrected p-Zonotopes of the state estimation error denoted by L i ∆ˆ x ,k and L i ∆¯ x ,k ,respectively, are given by L i ∆¯ x ,k = (cid:0) ∆ ¯ x k , G i ∆¯ x ,k , Σ i ∆¯ x ,k (cid:1) , L i ∆ˆ x ,k = (cid:0) ∆ ˆ x k , G i ∆ˆ x ,k , Σ i ∆ˆ x ,k (cid:1) , (8)where G i ∆¯ x and Σ i ∆¯ x denote the generator and covariance matrices of the corrected p-Zonotope of the state estimation error. Thecorrected p-Zonotope is estimated during the measurement update of set-valued DKF, which is explained later in Section 3.1.2.Similarly, G i ∆ˆ x and Σ i ∆ˆ x denote the generator and covariance matrices of the predicted p-Zonotope of the state estimation error.This predicted p-Zonotope is estimated during the time update of set-valued DKF, which is explained later in Section 3.1.1.Note that the above-defined covariance matrices Σ i ∆ˆ x and Σ i ∆¯ x represent the Gaussian characteristics of the p-Zonotopes andare different from ˆ P i,k and ¯ P i,k , defined earlier in Section 2.3.At the k th iteration, we utilize the translation set property in Eq. (2d) to represent the predicted and corrected p-Zonotopes ofthe global state vector L i ˆ x ,k and L i ¯ x ,k , respectively, as L i ˆ x ,k = x k + L i ∆ˆ x ,k and L i ∆¯ x ,k = x k + L i ∆¯ x ,k . Therefore, the predictedand corrected p-Zonotopes of the global state vector are given by Eq. (9). Note that translating the center of p-Zonotope asseen in Eq. (9) does not change the generator and covariance matrix. Therefore, the generator G and covariance matrix Σ of the predicted (or corrected) p-Zonotope of the state estimate x and state estimation error ∆ x are the same. From Eq. (9)and prior literature [16, 24], it is established that the center mean of predicted and corrected p-Zonotopes is equivalent to theirpoint-valued DKF estimates. L i ˆ x ,k = x k + L i ∆ˆ x ,k = (cid:0) ˆ x k , G i ∆ˆ x ,k , Σ i ∆ˆ x ,k (cid:1) L i ¯ x ,k = x k + L i ∆¯ x ,k = (cid:0) ¯ x k , G i ∆¯ x ,k , Σ i ∆¯ x ,k (cid:1) (9)8 .1.1 Predicted p-Zonotope of the state estimation error via time update In the time update of set-valued DKF, as shown in Fig. 5, we estimate the predicted p-Zonotope of the state estimation errorfor the next time iteration by utilizing the corresponding corrected p-Zonotope of the current iteration and the p-Zonotope ofprocess noise bounds. From the time update of point-valued DKF explained earlier in Eq. (5), we formulate the point-valuedstate estimation error ∆ ˆ x i,k as ∆ ˆ x i,k = ˆ x i,k − x k , = F ¯ x i,k − − ( F x i,k − + ν i,k ) , = F ∆ ¯ x i,k − − ν i,k . (a) Set-valued time update: Time bias (b) Set-valued time update: Time drift Figure 5: Visualization of the predicted p-Zonotope of the state estimation error after the time update of set-valued DKF. The2D p-Zonotope is sliced to be shown as two 1D p-Zonotopes, one each for (a) time bias; (b) time drift. The Minkowski sum ofthe linearly mapped predicted p-Zonotope at the k th iteration, indicated in orange, and the pre-computed process noise boundat the ( k − th iteration, indicated by gray, provide an estimate of the predicted p-Zonotope for the k th iteration, indicated bythe Parula colormap.Similar to the measurement update, given that ∆ ¯ x i,k − , ν ik ∀ j ∈ M i are independent random variables, we apply theproperties of p-Zonotopes discussed earlier in Eqs. (2a)-(2b) to estimate the predicted p-Zonotope as seen in Eqs. (10)-(11).For an example shown in Fig. 6, the predicted p-Zonotope obtained after the Minokwski sum in Eq. (10) is indicated by theParula colormap. L i ∆ˆ x ,k = F L i ∆¯ x ,k − (cid:77) L j ν,k , (10)such that G i ∆ˆ x ,k = (cid:20) F G i ∆¯ x ,k − , G i ν ,k (cid:21) , Σ i ∆ˆ x ,k = F Σ i ∆¯ x ,k − F (cid:62) + Σ i ν ,k . (11) In the measurement update step of set-valued DKF, as shown in Fig. 6, we utilize the p-Zonotope of measurement noise boundsand predicted p-Zonotope of the state estimation error to compute the corrected p-Zonotope. Unlike the conventional point-valued DKF that utilizes the optimal Kalman gain K j,k , we utilize an adaptive gain matrix ˜ K j,k to perform the measurementupdate of the set-valued DKF. This adaptive gain matrix is estimated during the spoofing attack mitigation module, which isexplained later in Section 3.2. From the measurement update of point-valued DKF explained earlier in Eq. (6), we formulate9he point-valued state estimation error ∆ ¯ x i,k as ∆ ¯ x i,k = ¯ x i,k − x k = ˆ x i,k + (cid:88) j ∈ M i ˜ K j,k (cid:16) y j,k − H j ˆ x i,k (cid:17) − x k = ∆ ˆ x i,k + (cid:88) j ∈ M i ˜ K j,k (cid:16) H j x k + ω j,k − H j ˆ x i,k (cid:17) ∆ ¯ x i,k = (cid:16) I − (cid:88) j ∈ M i ˜ K j,k H j (cid:17) ∆ ˆ x i,k + (cid:88) j ∈ M i ˜ K j,k ω j,k Given that ∆ ˆ x i,k , ω j,k ∀ j ∈ M i are independent random variables, we apply the properties of p-Zonotopes discussedearlier in Eqs. (2a)-(2b) to convert the point-valued representation of Eq. (6) to set-valued stochastic reachable sets as seen inEq. (12). For an example shown in Fig. 6, the corrected p-Zonotope obtained after the Minokwski sum in Eq. (12) is indicatedby a Parula colormap. L i ∆¯ x ,k = (cid:18) I − (cid:88) j ∈ Ω i ˜ K j,k H j (cid:19) L i ∆ˆ x ,k (cid:77) j ∈ M i ˜ K j,k L j ω (12)This is simplified further using Eq. (2c) to obtain Eq. (13). Intuitively, the estimated generator and covariance matrix, i.e., G i ∆¯ x ,k and Σ i ∆¯ x ,k provide bounds on the unknown state estimation error ∆ ¯ x i,k associated with the estimated point-valuedstate estimate ¯ x i,k . The corrected p-Zonotope L i ∆¯ x ,k of the state estimation error is analyzed to perform timing risk analysis inSection 3.3. G i ∆¯ x ,k = (cid:34)(cid:16) I − (cid:88) j ∈ M i ˜ K j,k H j (cid:17) G i ∆ˆ x ,k , ˜ K ,k G ω ,k , · · · , ˜ K (cid:12)(cid:12) M i (cid:12)(cid:12) ,k G (cid:12)(cid:12) M i (cid:12)(cid:12) ω ,k (cid:35) Σ i ∆¯ x ,k = (cid:16) I − (cid:88) j ∈ M i ˜ K j,k H j (cid:17) Σ i ∆ˆ x ,k (cid:16) I − (cid:88) j ∈ M i ˜ K j,k H j (cid:17) (cid:62) + (cid:88) j ∈ M i ˜ K j,k Σ j ω ˜ K (cid:62) j,k (13) (a) Set-valued measurement update: Time bias (b) Set-valued measurement update: Time drift Figure 6: Visualization of the corrected p-Zonotope of the state estimation error after the measurement update of set-valuedDKF. The 2D p-Zonotope is sliced to be shown as two 1D p-Zonotopes, one each for (a) time bias; (b) time drift. TheMinkowski sum of the linearly mapped predicted p-Zonotope at the k th iteration, indicated in red, and the pre-computedmeasurement noise bounds, indicated by gray, provide an estimate of the corrected p-Zonotope for the k th iteration, indicatedby the Parula colormap. Before measurement update at the i th receiver and k th iteration, we independently utilize the predicted p-Zonotope of the stateestimation error L i ∆ x ,k from Eq. (11) and the pre-computed measurement noise bound L iω,k to estimate the p-Zonotope of ex-pected measurement innovation. We apply the set properties of p-Zonotopes in Eqs. (2a)-(2c) on the point-valued measurement10nnovation (cid:15) i,k = z i,k − H i ˆ x i,k , which is defined earlier in Eq. (6). This p-Zonotope denoted by L i (cid:15) ,k is given by L i (cid:15) ,k = L i ω ,k (cid:77) H i L i ∆ˆ x ,k . (14)The p-Zonotope of expected measurement innovation L i (cid:15) ,k provides a probability value corresponding to each point-valuedmeasurement innovation (cid:15) i,k and this probability is denoted by α i,k . Based on this, for any i th receiver, we independentlycompute the attack status ˜ α i,k of the received GPS measurements by normalizing the obtained probability value α i,k with theprobability value corresponding to the center mean of the p-Zonotope L i (cid:15) ,k and subtracting it from one. The measurementattack status ˜ α i,k ∈ R lies between [0 , . As seen earlier in Eq. (7), the pre-computed process and measurement noise boundsare defined for non-spoofed conditions. Therefore, in a non-spoofed condition, the point-valued measurement innovation (cid:15) i,k and its p-Zonotope L i (cid:15) ,k of expected measurement innovation are in close agreement, and hence a low attack status ˜ α i,k ≈ isobtained. However, during spoofing, the point-valued measurement innovation does not comply with this estimated p-Zonotope,and therefore, a high value of ˜ α i,k ≈ is observed.Thereafter, each i th receiver in the network i ∈ M , broadcasts its received GPS measurement z i,k and the estimated attackstatus ˜ α i,k to all the receivers in its neighborhood set j ∈ M i − i . The set M i − i simply implies that the i th receiver need notbroadcast GPS measurements to itself. Similarly, each i th receiver also receives the GPS measurements z j,k and their associatedattack status ˜ α i,k from one or more receivers in its neighborhood set j ∈ M i − i . In an intuitive sense, the measurement attackstatus obtained from each j th receiver in the neighborhood set indicates its belief in the received GPS measurement z j,k .Thereafter, the measurement attack status is used to compute the adaptive gain matrix ˜ K j,k as seen in Eq. (15). Given that thegain matrix is a value between and , the Eq. (15) empirically represents a joint minimization of both spoofing probabilityand covariance of the point-valued state estimate. The formulated adaptive gain matrix ˜ K j,k is later utilized to perform theset-valued measurement update of SR-DKF, which was described earlier in Eq. (13). ˜ K ij,k = (cid:0) − ˜ α ij,k (cid:1) H (cid:62) j ¯ P i,k R − j,k ∀ j ∈ M i (15) We evaluate the probability of corrected p-Zonotope L i ∆¯ x,k of state estimation error derived in Eq. (13), to intersect an unsafeset, i.e., a set of states that violate the IEEE C37.118.1a-2014 standards. As shown in both Figs. 7(a) and 7(b), the unsafe setrepresents two gray strips, i.e., ∆ ¯ T i,k ∈ [26 . µs, ∞ ] and ∆ ¯ T i,k ∈ [ −∞ , − . µs ] . Intuitively, this matches the definition ofrisk described earlier in Section 1, which denotes the probability of the state estimation error to exceed a safety AL. (a) Over-approximation of corrected p-Zonotope (b) Top view of finite polytopes Figure 7: Visualization of the corrected p-Zonotope L i ∆¯ x,k of state estimation error and its intersection with the unsafe set B ;(a) shows an example that over-approximates the corrected p-Zonotope by four polytopes, which are denoted by D to D ;(b) shows the top-view of the finite polytopes and an example intersection region of D polytope with unsafe set B , which isindicated by red cross marks.To perform the timing risk analysis, we over-approximate the corrected p-Zonotope L i ∆¯ x,k by a finite series of poly-topes [24]. An example of over-approximation of p-Zonotope by four polytopes D to D is shown in Fig. 7(a). There-after, we compute the intersection area of each l th polytope D l with the unsafe set B . Figure 7(b) shows the top-view ofover-approximated series of polytopes and an example intersection region of D polytope with unsafe set B .11he timing risk, seen in Eq. (16), is defined as the sum across polytopes D l of their intersection area, which is denotedby A ( · ) , with the unsafe set B times the max probability of D l th polytope, which is denoted by p max,l . The max probability of D l th polytope is at its top face and equals the probability of corrected p-Zonotope at that level. For computational tractability,we represent the p-Zonotope by a finite number of polytopes and a γ -confidence set. The γ -confidence set [24] thresholds thetail of a p-Zonotope by a pre-defined parameter γ , such that P [ − γ < N (0 , < γ ] = erf ( γ √ , where erf denotes the errorfunction [30]. From [24], the probability that the state estimation error lies outside the γ -confidence set of corrected p-Zonotopeis − erf (cid:16) γ (cid:112) (2) (cid:17) n . r i,k = 1 − erf (cid:16) γ √ (cid:17) n + (cid:88) l A ( D l ∩ B ) × p max,l (16) Figure 8 summarizes the implementation steps of the proposed SR-DKF described in Sections 3.1-3.3. We assume that the i threceiver ∀ i ∈ M is authentic during initialization. Therefore, at iteration k = 0 , we define the predicted p-Zonotope L i ∆¯ x , − of the state estimation error to be zero-mean center, with non-zero generator and over-bounding covariance matrices.Figure 8: Implementation steps of the proposed SR-DKF algorithm. We validate the robustness of the proposed SR-DKF algorithm to not only mitigate the effect of spoofing but also estimatethe secure GPS timing and its associated timing risk across all the receivers in the network. We perform set operations ofstochastic reachability and transform across various set representations, such as polytopes, zonotopes, p-Zonotopes, etc., usinga publicly-available MATLAB toolbox, known as COntinuous Reachability Analyzer (CORA) [31].As shown in Fig. 9, we design a simulated network of seven GPS receivers, (cid:12)(cid:12) M (cid:12)(cid:12) = 7 , that are spatially distributed acrossthe US. The network of GPS receivers indexed from Rx : 1 to Rx : 7 are located in 1) Stanford, CA; 2) Atlanta, GA; 3) Urbana,IL; 4) Boulder, CO; 5) Austin, TX; 6) Boston, MA; and 7) Auburn, AL. We utilized a C++ language-based software-definedGPS simulator known as GPS-SIM-SDR [32, 33] to simulate the raw GPS signals received at each location in the network.Later, we post-processed the simulated GPS signals using a MATLAB-based software-defined radio known as SoftGNSS [34].We define the simulated errors observed in the system during authentic conditions to have the following characteristics:1. For the process noise in GPS time (the first global state T i,k ), the mean and covariance of the time-varying Gaussiandistribution lies between [ − . µ s, . µ s] and [ µ s, µ s], respectively; Similarly, for the process noise in its drift rate(the second global state ˙ T i,k ), the time-varying mean and covariance lies between [ − . ns/s, . ns/s] and [ ns/s, ns/s],respectively; 12. For the measurement noise in GPS pseudoranges, the mean and covariance of the time-varying Gaussian distributionlies between [ − µ s, µ s] and [ µ s, µ s], respectively; Similarly, for the measurement noise in GPS Doppler, thetime-varying mean and covariance lies between [ − . ns/s, . ns/s] and [ ns/s, ns/s], respectively. We consider all thereceivers in this simulated setup to have the same bounds.3. For the initial state error bounds in GPS time (the first global state T i,k ), the mean and covariance of the time-varyingGaussian distribution lies between [ − . µ s, . µ s] and [ µ s, µ s], respectively; Similarly, for the initial state errorbounds in the second global state ˙ T i,k , the time-varying mean and covariance lies between [ − . ns/s, . ns/s] and[ ns/s, ns/s], respectively;As discussed earlier in Fig. 8, we initialize our SR-DKF algorithm by pre-defining the following p-Zonotopes ∀ i ∈ M :process noise L i ν , measurement noise L i ω , and initial state estimation error L i ∆¯ x , − . We formulate the p-Zonotopes from theabove-listed errors bounds in the same manner as the example discussed earlier in Section 2.1.Figure 9: Geographical map showcasing a simulated network of seven spatially distributed GPS receivers across the US. In the first set of experiments, we consider sparse connectivity among the seven receivers in the network. Figure 10(a) showcasesthe connectivity matrix where a green tick indicates the presence of a bidirectional communication link while a red crossindicates otherwise. Note that each receiver is connected to itself by default. We induce a simulated coordinated signal-levelspoofing attack in the simulated GPS measurements received at two GPS receivers in the network, i.e., Rx : 1 , which is locatedin Stanford, CA and Rx : 5 , which is located in Austin, TX. Between the time duration k = 40 − s, the first spoofingattack manipulates the GPS time at Rx : 5 to deviate at a rate of ns/s, and the second spoofing attack, which occurs between k = 800 − s, deviates the GPS time of Rx : 1 at a rate of ns/s.In Fig. 10(b), we show the performance of the point-valued state estimation via an adaptive DKF framework, where even thereceivers that are not directly spoofed, e.g., Rx : 3 , exhibit high state estimation errors. Also, in Fig. 10(c), we show the stateestimation accuracy achieved via a point-valued single-receiver-based adaptive KF. In both these point-valued approaches, themeasurement noise covariance R i,k is adaptively estimated via Eq. (6) by setting the forgetting factor ψ = 0 . . Based on themaximum error values listed in Table 1, we observe that both the approaches violate the IEEE C37.118.1a-2014 standards. InFig. 10(d), we validate that our proposed SR-DKF algorithm provides a secure estimate of GPS timing with maximum error of . µ s across all receivers, and thereby complies with the IEEE standards at all times. The quantitative maximum error valuesof the receiver network for the above-described state estimation approaches are listed in Table 1.13 a) Sparse connectivity across the network (b) Point-valued state estimation via adaptive DKF(c) Point-valued adaptive KF (d) Set-valued state estimation via SR-DKF Figure 10: Performance comparison of the three state estimation techniques; a) sparse connectivity across the network ofreceivers; b) set-valued state estimation using the proposed SR-DKF algorthm; c) point-valued adaptive DKF; d) point-valuedadaptive KF for a single receiver; In (b)-(d) the vertical dotted black lines indicate the start and stop times of spoofing attackswhile the horizontal magenta dotted lines indicate the threshold set by the IEEE C37.118.1a-2014 standards.Table 1: Comparison among the maximum error values of the global state vector estimated via three approaches: the proposedSR-DKF algorithm, point-valued adaptive DKF, and point-valued single-receiver-based adaptive KF.State estimation approach Rx : 1 Rx : 2 Rx : 3 Rx : 4 Rx : 5 Rx : 6 Rx : 7 Set-valuedSR-DKF 8.27 7.42 6.81 6.79 8.43 7.10 7.40
Point-valuedadaptive DKF 35.88 6.47 15.24 6.70 74.68 6.26 6.36 T ( µ s) Point-valuedadaptive KF 257.36 7.36 9.48 6.43 118.90 6.46 6.48 Set-valuedSR-DKF 22.05 19.12 20.54 17.12 33.23 16.98 18.01
Point-valuedadaptive DKF 80.68 23.60 39.94 19.95 125.23 18.01 20.88 ˙ T (ns/s) Point-valuedadaptive KF 415.39 17.90 22.75 18.71 122.94 16.21 17.04In Fig. 11, we plot the attack status ˜ α i,k of the GPS measurements z i,k at each receiver in the network. In accordancewith Section 3.2, we demonstrate that the proposed SR-DKF algorithm successfully detects and mitigates spoofing in the GPSmeasurements of Rx : 1 and while accurately categorizing the other receivers, i.e., Rx : 2 , , , , to be non-spoofed withattack status ˜ α i,k ≈ at all times. We observe that the spoofing attack with higher drift rate is mitigated quicker, i.e., ns/s14n Rx : 1 , as compared to the one with lower drift rate, i.e., ns/s in Rx : 5 .In Fig. 12, we showcase the timing risk associated with the estimated GPS timing, i.e., the probability of GPS time to violatethe IEEE C37.118.1a-2014 standards. This is computed by analyzing the intersection of the corrected p-Zonotope L i ∆¯ x ,k withthe unsafe set B , as discussed in Section 3.3. An increase or jump in measure of timing risk is based on the adaptive gainmatrix ˜ K j,k ∀ j ∈ M i that in turn depends on the estimated attack status ˜ α j,k of the receiver. Therefore, an increase in thetiming risk for Rx : 5 till t = 1040 s is consistent with the decrease in measurement reliability, which is evident from Fig 11.By comparing Fig. 10(a) and Fig. 12, we observe a direct correlation between the value of timing risk to the number of reliablemeasurements available for the calculation of corrected p-Zonotope. For instance, since none of the four neighbors of Rx : 4 are spoofed, the estimated timing risk is consistently low, i.e., ≈ − , for the entire experiment duration. In contrast, Rx : 3 has significantly high timing risk ( ≈ . because two out of its three neighbors are spoofed. Similarly, when only one out ofthree neighbors at Rx : 1 is spoofed, the order of timing risk magnitude is − whereas when one out of two neighbors arespoofed at Rx : 5 , the timing risk is a higher order magnitude of . . Therefore, our SR-DKF estimates a robust measure of thetiming risk across all receivers in the network.Figure 11: Attack status of each receiver in the network that is estimated via the proposed SR-DKF algorithm. This is computedby evaluating the point-valued measurement innovation against the estimated p-Zonotope of expected measurement innovation,as seen in Eq. (14). Our SR-DKF adaptively mitigates spoofing in Rx : 1 , and successfully estimates other receivers to benon-spoofed.Figure 12: Robust timing risk analysis for each receiver in the network via the proposed SR-DKF algorithm. This is computedby evaluating the intersection of the corrected p-Zonotope of the state estimation error with the unsafe set as seen in Eq. (16).Our SR-DKF successfully quantifies the risk of violating the IEEE C37.118.1a-2014 standards.15 .1.1 Robustness analysis of the proposed SR-DKF algorithm In the second set of experiments, we demonstrate the robustness of the proposed SR-DKF algorithm as the number of receiversin the neighborhood set | M | is varied. Considering a total experiment duration of s and a network of seven receivers seenin Fig. 9, we induced a simulated meaconing attack between k = 10 − s in the GPS measurements that is received at Rx : 1 located at Stanford, CA.In Fig. 13, we analyzed the variation in the timing risk associated with the estimated GPS time of Rx : 1 as the numberof receivers in the network M is increased. For each size of the receiver network that ranges from | M | = 2 to | M | = 7 , weconsider the presence of a fully-connected bidirectional communication network and perform Monte-Carlo runs for fourmeaconing magnitudes of µ s, µ s, µ s, and µ s. We successfully demonstrated that as the number of receiversinteracting with the Rx : 1 , i.e., | M | , increases, the associated timing risk decreases. Particularly, for | M | > , the decreasein timing risk is on the order of ≤ − and hence, considered insignificant. Similarly, the effect of meaconing magnitude onthe timing risk decreases as the number of interacting receivers associated with Rx : 1 increases. Therefore, we validated theproposed SR-DKF algorithm that requires only a handful of distributed GPS receivers to achieve robust performance.Figure 13: Variation in the timing risk associated with GPS timing of Rx : 1 with the number of receivers in a fully-connectednetwork. The colors red, blue, cyan and magenta indicate a meaconing attack of µ s, µ s, µ s, and µ s, respectively.We observe that the variation in timing risk decreases with an increase in the number of interacting receivers with Rx : 1 andthe magnitude of meaconing attack. In summary, we developed a Stochastic Reachability-based Distributed Kalman Filter (SR-DKF) to perform set-valued stateestimation using a network of GPS receivers that estimate not only the point-valued mean but also the stochastic reachable set ofstate (i.e., GPS time and its drift rate). By parameterizing the stochastic reachable set via probabilistic zonotopes (p-Zonotopes),we derived the time and measurement update steps of the set-valued DKF to estimate the predicted and corrected p-Zonotopeof timing error, respectively. To compute secure GPS time for a network of receivers, we designed a two-tiered approach: first,we performed spoofing mitigation at measurement-level via deviation of measurement innovation from its expected stochasticset (represented by the p-Zonotope). Second, we performed timing risk analysis at state-level via intersection probability ofcorrected p-Zonotope with an unsafe set, i.e., a set that violates the IEEE C37.118.1a-2014 standards. To our knowledge, noprior literature exists on secure state estimation via stochastic reachability nor on timing risk quantification during spoofing.Under a simulated coordinated signal-level spoofing that affects two among a network of seven sparsely-connected GPSreceivers, the proposed SR-DKF algorithm demonstrated a higher timing accuracy, successfully mitigated the spoofing attacks,and estimated a robust measure of timing risk. While our algorithm showed a maximum timing error of only < . µ s, both theconventional DKF and single-receiver-based adaptive KF showed high estimation errors that violate the IEEE C37.118.1a-2014threshold of . µ s. Furthermore, by varying the number of receivers and attack magnitudes, we performed extensive MonteCarlo runs to validate the sensitivity of estimated timing risk. 16 CKNOWLEDGEMENTS
We would like to thank Ashwin Kanhere, Siddharth Tanwar, and Tara Yasmin Mina for their help in revising the paper andpresentation. We would also like to thank the members of NAVLab research at Stanford University for their thoughtful inputon this work and paper.This material is based upon work supported by the Department of Energy under Award Number DE-OE0000780.This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the UnitedStates Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes anylegal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or processdisclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercialproduct, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply itsendorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions ofauthors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. R EFERENCES [1] G. Duan, Y. Yan, X. Xie, H. Tao, D. Yang, L. Wang, K. Zhao, and J. Li, “Development status quo and tendency of widearea phasor measuring technology,”
Automation of Electric Power Systems , vol. 39, no. 1, pp. 73–80, 2015.[2] P. Overholt, D. Ortiz, and A. Silverstein, “Synchrophasor technology and the doe: Exciting opportunities lie ahead indevelopment and deployment,”
IEEE Power and Energy Magazine , vol. 13, no. 5, pp. 14–17, 2015.[3] J. Taft, “Assessment of existing synchrophasor networks,”
Pacific northwest national laboratory (PNNL) technical report.https://gridarchitecture. pnnl. gov/media/whitepapers/Synchrophasor net assessment final. pdf. Accessed April , 2018.[4] P. T. Myrda and K. Koellner, “NASPInet - the internet for synchrophasors,” in , 2010, pp. 1–6.[5] L. Heng, D. B. Work, and G. X. Gao, “GPS signal authentication from cooperative peers,”
IEEE Transactions on Intelli-gent Transportation Systems , vol. 16, no. 4, pp. 1794–1805, 2014.[6] P. Misra and P. Enge, “Global positioning system: signals, measurements and performance second edition,”
Global Posi-tioning System: Signals, Measurements And Performance Second Editions , vol. 206, 2006.[7] D. P. Shepard, T. E. Humphreys, and A. A. Fansler, “Evaluation of the vulnerability of phasor measurement units to GPSspoofing attacks,”
International Journal of Critical Infrastructure Protection , vol. 5, no. 3-4, pp. 146–153, 2012.[8] S. Bhamidipati and G. X. Gao, “GPS multi-receiver joint direct time estimation and spoofer localization,”
IEEE Transac-tions on Aerospace and Electronic Systems , vol. 55, no. 4, pp. 1907–1919, 2018.[9] K. E. Martin, “Synchrophasor measurements under the IEEE standard C37. 118.1-2011 with amendment C37. 118.1 a,”
IEEE Transactions on Power Delivery
Proceedings of the 2014 ACM conference on Security and privacy in wireless & mobile networks , 2014, pp.99–104.[12] A. Khalajmehrabadi, N. Gatsis, D. Akopian, and A. F. Taha, “Real-time rejection and mitigation of time synchronizationattacks on the global positioning system,”
IEEE Transactions on Industrial Electronics , vol. 65, no. 8, pp. 6425–6435,2018.[13] S. Bhamidipati, T. Y. Mina, and G. X. Gao, “GPS time authentication against spoofing via a network of receivers for powersystems,” in . IEEE, 2018, pp. 1485–1491.1714] C. G¨unther, “A survey of spoofing and counter-measures,”
NAVIGATION: Journal of The Institute of Navigation , vol. 61,no. 3, pp. 159–177, 2014.[15] D. Schmidt, K. Radke, S. Camtepe, E. Foo, and M. Ren, “A survey and analysis of the GNSS spoofing threat and coun-termeasures,”
ACM Computing Surveys (CSUR) , vol. 48, no. 4, pp. 1–31, 2016.[16] D. Shi, T. Chen, and L. Shi, “On set-valued Kalman filtering and its application to event-based state estimation,”
IEEETransactions on Automatic Control , vol. 60, no. 5, pp. 1275–1290, 2014.[17] V. Shiryaev and E. Podivilova, “Set-valued estimation of linear dynamical system state when disturbance is decomposedas a system of functions,”
Procedia Engineering , vol. 129, pp. 252–258, 2015.[18] P. L. Combettes, “The foundations of set theoretic estimation,”
Proceedings of the IEEE , vol. 81, no. 2, pp. 182–208,1993.[19] M. Prandini, J. Hu, C. Cassandras, and J. Lygeros, “Stochastic reachability: Theory and numerical approximation,”
Stochastic hybrid systems, Automation and Control Engineering Series , vol. 24, pp. 107–138, 2006.[20] S. Bhamidipati and G. X. Gao, “Multiple GPS fault detection and isolation using a graph-SLAM framework,” in . Institute ofNavigation, 2018, pp. 2672–2681.[21] ——, “GPS spoofing mitigation and timing risk analysis in networked PMUs via stochastic reachability,” in . Institute ofNavigation (Accepted), 2020.[22] S. P. Talebi and S. Werner, “Distributed Kalman filtering: Consensus, diffusion, and mixed,” in . IEEE, 2018, pp. 1126–1132.[23] A. Alanwar, H. Said, and M. Althoff, “Distributed secure state estimation using diffusion Kalman filters and reachabilityanalysis,” in . IEEE, 2019, pp. 4133–4139.[24] M. Althoff, O. Stursberg, and M. Buss, “Safety assessment for stochastic linear systems using enclosing hulls of proba-bility density functions,” in . IEEE, 2009, pp. 625–630.[25] J. R. Nu˜nez, C. R. Anderton, and R. S. Renslow, “Optimizing colormaps with consideration for color vision deficiency toenable accurate interpretation of scientific data,”
PloS one , vol. 13, no. 7, p. e0199239, 2018.[26] M. Althoff and J. M. Dolan, “Online verification of automated road vehicles using reachability analysis,”
IEEE Transac-tions on Robotics , vol. 30, no. 4, pp. 903–918, 2014.[27] E. W. Grafarend,
Linear and nonlinear models: fixed effects, random effects, and mixed models . de Gruyter, 2006.[28] L. Heng, G. X. Gao, T. Walter, and P. Enge, “GPS signal-in-space anomalies in the last decade: Data mining of400,000,000 GPS navigation messages,” in , 2010, pp. 3115–3122.[29] S. T. Watt, S. Achanta, H. Abubakari, E. Sagen, Z. Korkmaz, and H. Ahmed, “Understanding and applying precision timeprotocol,” in . IEEE, 2015, pp. 1–7.[30] K. B. Oldham, J. C. Myland, and J. Spanier, “The error function erf (x) and its complement erfc (x),” in
An Atlas ofFunctions . Springer, 2008, pp. 405–415.[31] M. Althoff, “CORA 2016 manual,”
TU Munich , vol. 85748, 2016.[32] T. Ebinuma, “GPS-SDR-SIM,” [Online] Available: https://github.com/osqzss/gps-sdr-sim, (Accessed as of June, 2020).[33] S. Bhamidipati, K. J. Kim, H. Sun, and P. V. Orlik, “GPS spoofing detection and mitigation in PMUs using distributedmultiple directional antennas,” in