GPS Spoofing Mitigation and Timing Risk Analysis in Networked Phasor Measurement Units via Stochastic Reachability

Sriramya Bhamidipati; Grace Gao

doi:10.33012/navi.574

Abstract

To address phasor measurement unit (PMU) vulnerability to spoofing, we propose the use of a set-valued state estimation technique known as stochastic reachability (SR)-based distributed Kalman filter (DKF) that computes secure global positioning system (GPS) timing across a network of receivers. Utilizing SR, we estimate not only GPS time but also its stochastic reachable set, which is parameterized by probabilistic zonotope (p-Zonotope). While requiring known measurement error bounds in only non-spoofed conditions, we designed a two-tiered approach. We first performed measurement-level spoofing mitigation via deviation of a measurement innovation from its expected p-Zonotope. We then performed state-level timing risk analysis via a determination of the intersection probability of the estimated p-Zonotope with an unsafe set that violates IEEE C37.118.1a-2014 standards. Finally, we validated our SR-DKF algorithm by subjecting it to a simulated receiver network to coordinate signal-level spoofing. We demonstrate improved timing accuracy and successful spoofing mitigation via the use of our SR-DKF algorithm. We also validated the robustness of the estimated timing risk as the number of receivers were varied.

Keywords

1 INTRODUCTION

Modern power systems analyze real-time grid stability using a widely-dispersed network of advanced devices known as phasor measurement units (PMUs) (Al-Hammouri et al., 2012). As per the American Recovery and Reinvestment Act of 2009, more than 2500 PMUs were deployed across North America as of March 2015 (Tholomier et al., 2009). Figure 1 shows a map of the dense network of these spatially-distributed PMUs (Dagle, 2018). The North American synchrophasor initiative network (NASPInet), which is a collaborative effort by the United States Department of Energy and the North American Electric Reliability Corporation, provides a secure and decentralized infrastructure for communication among the PMUs (Myrda & Koellner, 2010). Each PMU is equipped with a Global Positioning System (GPS) receiver that precisely time-stamps the voltage and current phasor values at critical power substations. A network of geographically distributed GPS receivers facilitates robust time authentication by internal cross-checking, thereby eliminating the need for an external reference receiver (Heng et al., 2015). The networked approach in power systems overlays the existing grid infrastructure and can therefore be implemented without the need for additional hardware.

FIGURE 1

A network of more than 2500 PMUs deployed across North America as of March 2015 (Dagle, 2018)

PMUs are susceptible to GPS spoofing because of the low signal power and unencrypted structure of civilian GPS/L1 signals (Hofmann-Wellenhof et al., 1992). As validated in prior literature (Shepard et al., 2012), this vulnerability can be exploited by attackers who can mislead the GPS timing at one or more critical substations and cause cascading system failure. The various spoofing attacks that have the potential to disrupt GPS timing vary from simple meaconing to more sophisticated signal-level spoofing and coordinated attacks. A spoofer executing meaconing (or time jump) broadcasts earlier recorded authentic GPS signals which are then used to induce a constant lag in the victim receiver’s computed time compared to the true time. In signal-level spoofing (or time walk), the spoofer broadcasts the simulated GPS signals and will gradually steer the computed time of victim receiver away from its true time. Coordinated attack occurs when multiple spoofers manipulate the GPS signals received at one or more victim receivers. We refer readers to our earlier work (Bhamidipati & Gao, 2019) for a detailed overview of these types of attacks.

Secure GPS timing is of paramount importance for reliable PMU operations. This requires ensuring that the timing accuracy satisfies the IEEE C37.118.1a-2014 standard “Synchrophasors for Power Systems” (Martin, 2015). The IEEE C37.118.1a-2014 standard defines a metric called total vector error (TVE) which is the square root of the difference squared between the real and imaginary parts of the theoretical phasor and the estimated phasor to the ratio of magnitude of the theoretical phasor and presented as a percentage. An error in estimated timing presents as an error in the phase angle of voltage and current phasors. Furthermore, without magnitude errors, a phase angle error of 0.573° ≈ 0.01 rad corresponds to a 1% TVE, which is the maximum allowable TVE by the IEEE C37.118.1a-2014 Standard. To ensure that the phase angle error is less than 0.01 rad, we require the error in GPS timing, ΔT, to satisfy 2πf_sampΔT < 0.01 rad, where f_samp is the PMU’s data logging frequency. For a PMU operating at f_samp = 60 Hz, which is a commonly used frequency for phasor data logging (Martin, 2015), less than 1% TVE corresponds to a GPS timing error |ΔT| ≤ 26.5 μs.

To address concerns regarding GPS vulnerability against spoofing in power systems, accounting for timing accuracy as a standalone metric will not be sufficient. It is equally important to quantify the risk associated with the estimated timing. While there are different ways (Georgescu et al., 2019; Ochieng et al., 2003) to quantify risk, for this work, we define risk as the probability of an estimation error that exceeds a safety alert limit (AL). For instance, claiming that the estimated timing is 1 μs accurate does not suffice, unless this is accompanied by a certain risk, say 10⁻³. In this work, we ensure secure GPS timing in power systems by performing not only spoofing attack mitigation but also timing risk analysis. However, there are no defined standards for acceptable risk in power systems at this time. Therefore, we set a conservative AL = ±1 μs, as compared to the aforementioned IEEE C37.118.1a-2014 standard. This is done to ensure that the estimated GPS time complies with the commonly desired timing accuracy (Heng et al., 2015) and never exceeds the IEEE C37.118.1a-2014 threshold of ±26.5 μs without raising an alert.

The state-of-the-art prior work (Bhamidipati et al., 2018; Khalajmehrabadi et al., 2018; Risbud et al., 2018) on networked anti-spoofing methods falls under the category of point-valued state estimation approaches. A point-valued state estimate denotes the expected value of the state vector for a given set of measurements. The existing literature focused on of point-valued anti-spoofing methods (Bhamidipati & Gao, 2019; Günther, 2014; Schmidt et al., 2016) includes studies that discuss signal quality monitoring, angle-of-arrival, and time-difference-of-arrival-based analyses of the incoming satellite signals, validation of the encryption codes, 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 (Shi et al., 2015). An illustration of the difference between point-valued and set-valued propagation is shown in Figure 2(a).

FIGURE 2

Set-valued approach and stochastic reachability. Shown in (a) is an illustration of the difference between point-valued and set-valued approaches (non-stochastic). Unlike the point-valued techniques in which only the point-valued state is propagated across time, the set-valued approach propagates a set of state estimates, thereby providing bounds on the estimation errors. The image in (b) provides an intuitive understanding of stochastic reachability that computes a set of stochastic reachable states given an initial stochastic set and known error bounds.

There is much prior available literature that focuses on set-valued theory (Shi et al., 2015; Shiryaev & Podivilova, 2015). Our current work was inspired by a set-valued formulation designed for stochastic processes known as stochastic reachability (Bujorianu, 2012). While stochastic reachability is widely used for path planning and collision avoidance in robotics applications, it has not been applied in the field of power systems. As shown in Figure 2(b), stochastic reachability computes a set of stochastic reachable states given an initial stochastic set of states and known stochastic sets of system disturbances and measurement errors. In this context, each state included in the stochastic reachable set is associated with a measure that indicates its maximum probability of occurrence. The use of stochastic reachability ensures secure state estimation while attempting to account for difficult-to-model and explicitly neglected process or measurement errors in the receiver clock dynamics, communication latency, and spoofed GPS measurements (i.e., when the associated error distribution is not known). For brevity, these errors will be referred to as unmodeled biases.

In this work, we propose a set-valued state estimation technique that can be applied to a network of GPS receivers, termed Stochastic Reachability-Based Distributed Kalman Filter (SR-DKF). While the stochastic reachability formulation requires known error bounds, the attack errors introduced by GPS spoofing are unbounded and time-varying in nature. To account for this, we consider the error bounds in GPS measurements to be known factors only in non-spoofed (i.e., authentic) conditions. We pre-compute these bounds via offline empirical analysis of historical data (Bhamidipati & Gao, 2018). We parameterize the stochastic reachable set via an enclosing probabilistic hull known as probabilistic Zonotope (p-Zonotope) (Althoff et al., 2009) that can efficiently model stochastic biases. This work is based on our recent ION GNSS+2020 conference paper (Bhamidipati & Gao, 2020). Note that we modified the mathematical notations compared to those presented in our prior work (Bhamidipati & Gao, 2020) to ensure consistency with the standard notations utilized in estimation theory (Crassidis & Junkins, 2011; Simon, 2007).

In a conventional point-valued DKF approach (Talebi & Werner, 2018), each receiver compares the deviation of the received measurements against those that are expected and thereafter fuses the measurements obtained from different receivers to compute a point-valued state estimate. However, this approach is designed for use under nominal fault-free conditions with the initial states, process, and measurement noise modeled to follow a Gaussian distribution. Therefore, this approach cannot account directly for unmodeled biases in the system. Other point-valued DKF approaches (Li et al., 2013; Najafi et al., 2019) account for non-zero biases by adaptively estimating the process and measurement noise covariances. However, approximating unmodeled biases using a zero-mean Gaussian distribution can become restrictive, thereby limiting the practical applicability of these approaches. Of note, these point-valued techniques must be fused with independent modules, i.e., fault mitigation and integrity monitoring, in order to detect anomalies and provide a measure of timing risk. Our proposed SR-DKF describes an alternative approach to previously-identified methods that focus on point-valued anti-spoofing. In contrast to the point-valued techniques, our proposed SR-DKF algorithm performs secure state estimations of GPS time and spoofing attack mitigation by propagating these unmodeled biases (i.e., those with known bounds but unknown distribution) in the set-valued domain. Therefore, we estimate not only the point-valued mean but also the stochastic reachable set of GPS times that are later utilized to compute the associated risk of violating the pre-defined safety AL. Essentially, our proposed SR-DKF algorithm combines the three important steps, namely, secure state estimation, fault mitigation, and timing risk analysis into a single integrated pipeline, thus ensuring tolerance against anomalies. Essentially, our proposed SR-DKF algorithm provides a sequential estimate of stochastic timing output that remains valid even in the presence of anomalies, wherein the validity is quantified based on the estimated measure of timing risk. To the best of our knowledge, there is no prior literature that addresses methods used to perform secure state estimation using stochastic reachability nor are there any that focus on quantifying the timing risk associated with the estimated GPS time.

1.1 Design Aspects of the Proposed SR-DKF Algorithm

Here we present the three main design aspects of the proposed SR-DKF algorithm, including i) a spatially-distributed network of static GPS receivers, ii) a distributed processing framework, and iii) the set representation via p-Zonotopes, which are described as follows:

1.1.1 A spatially-distributed network of static GPS receivers

We utilize the existing grid infrastructure shown earlier in Figure 1 to consider a small group of spatially-distributed GPS receivers (approximately four to eight). This framework provides two advantages. First, by considering a spatially-distributed network, we increase the measurement redundancy while minimizing the probability of multiple receivers being simultaneously spoofed. Second, given that the power substations are static, we can utilize the known receiver positions to aid the proposed SR-DKF algorithm.

1.1.2 A distributed processing framework

We opted for distributed processing to minimize the reaction time to attacks and also to localize the extent of failure. This is justified because the GPS reference time for synchronizing the PMUs is no longer reliable during spoofing. Each receiver in our distributed network asynchronously broadcasts the latest GPS measurement data to its neighbors, thereby reducing communication latency.

We consider a combination of an independent clock model at each power substation and communication links between various power substations to provide short-term time synchronization within the network. By contrast, GPS receiver algorithms (in this work, our proposed SR-DKF algorithm) play a key role in maintaining a bound on both the synchronization errors between any two receivers and the absolute time at each receiver in the network. This indicates a standard implementation of PMUs/power grids since the development of newer clocks, for example, the chip-scale atomic clock (CSAC) and communication protocols such as network-time protocol (NTP) and precise-time protocol (PTP) that can provide up to tens of microsecond-level accuracy over shorter time spans (Allnutt et al., 2017; Wang et al., 2008; Zhan et al., 2017).

1.1.3 A set representation via p-Zonotopes

Among all available set representations (Althoff & Dolan, 2014) such as zonotope, ellipsoid, polytope, and others, our choice of p-Zonotope was justified because it can efficiently enclose the stochastic biases in GPS measurements. A p-Zonotope can enclose a variety of distributions such as mixture models, distributions with uncertain/time-dependent mean, and non-Gaussian/unknown distributions, among others. A one-dimensional (1D) example of a p-Zonotope that encloses multiple Gaussian distributions with varying mean and covariance is illustrated in Figure 3(a). The probabilistic nature of p-Zonotopes blends well with the Gaussian-centric framework of DKF.

FIGURE 3

Illustration of p-Zonotopes (a) shows an example of a 1D p-Zonotope as indicated by the Parula colormap (Nuñez et al., 2018). The range of colors (from blue to yellow) denotes the lowest to highest values of probability represented by this p-Zonotope. This 1D p-Zonotope encompasses multiple Gaussian distributions, a subset of which are shown in gray. (b) is an example of a 2D p-Zonotope as indicated by the Parula colormap. This 2D p-Zonotope can be used to represent the state vector, which is explained in Section 2.

Unlike the conventional techniques that propagate point-valued state estimates, the proposed SR-DKF algorithm utilizes pre-computed error bounds on the receiver clock dynamics and measurements to propagate the p-Zonotope of the state estimate recursively. An example of a two-dimensional (2D) p-Zonotope is shown in Figure 3(b). In this work, the 2D p-Zonotope is used to represent the state estimate that is introduced later in Section 2.

1.2 Contributions

The main contributions of this paper are as follows:

Utilizing an efficient set representation known as p-Zonotope (Althoff et al., 2009) to parameterize the stochastic reachable set and account for stochastic biases, we develop a set-valued DKF framework to compute secure GPS time for a network of receivers.
Considering offline-estimated measurement error bounds in non-spoofed (i.e., authentic) conditions, we design a two-tier approach that performs:
1. Measurement-level spoofing attack mitigation by adaptive evaluation of the deviation of point-valued GPS measurements from the set-valued domain of measurement innovation. Therefore, the proposed SR-DKF algorithm is resilient to a variety of GPS attacks, including meaconing, signal-level spoofing, and coordinated attacks, among others.
2. State-level timing risk analysis by estimating the intersection probability of the estimated set of stochastic reachable states (represented by p-Zonotope) with an unsafe set, i.e., a set that violates the pre-defined safety AL (Martin, 2015).
For a simulated setup, we demonstrate the successful mitigation of both simple meaconing and sophisticated signal-level spoofing attacks. During a simulated case of coordinated attack that affects multiple receivers, we show an improved state estimation accuracy via the proposed SR-DKF algorithm compared to conventional point-valued and adaptive DKF. We also demonstrate the robustness of the proposed SR-DKF algorithm by analyzing the variation in timing risk with respect to the number of receivers and the magnitude of spoofing attacks.

The remaining sections of the paper are organized as follows: Section 2 describes the preliminaries of p-Zonotopes; Section 3 outlines the system model; Section 4 describes the proposed SR-DKF algorithm; Section 5 experimentally validates the improved state estimation accuracy of GPS time computed via the proposed SR-DKF algorithm as compared to that estimated by point-valued DKF techniques, and also computes a robust measure of timing risk associated with the estimated GPS time; and Section 6 concludes the paper.

2 PRELIMINARIES OF P-ZONOTOPES

A p-Zonotope , as seen in Equation (1) and shown earlier in Figure 3(a), is characterized by three parameters (Althoff et al., 2009), namely, a) the mean of its center, denoted by , b) the uncertainty in the center, represented by the generator matrix , and c) the covariance of Gaussian tails, denoted by .

Intuitively, a p-Zonotope with a certain center mean, i.e., zero center uncertainty G = 0, represents a Gaussian distribution that overbounds multiple probability density functions (PDFs) with the same center mean c. This is represented by in Equation (2).

where {·} denotes the set representation, denotes pairwise independent Gaussian distributed random variables, and denotes the associated generator vectors of , such that with [·;·] as the horizontal concatenation operator. Note that, the covariance ∑ of a p-Zonotope defined earlier in Equation (1) is equivalent to .

Similarly, a p-Zonotope with zero covariance Σ = 0 simplifies to a zonotope (Althoff & Dolan, 2014) that encompasses all the possible values of the center, such that the generator matrix of Z is G = [g¹, ⋯, g^e]. This case is mathematically represented by Z in Equation (3).

where denotes independent numbers. Note that the parameter G defined in Equations (1) and (3) is not the same as defined in Equation (2).

The p-Zonotope is a combination of both sets, i.e., , where is a set-addition operator. Alternatively, , where represents the PDF of distributions represented by and denotes the expectation operator. From Equations (1)–(3), note that depends on both G and . Additionally, unlike a PDF, the area enclosed by a p-Zonotope is not equal to one. More details regarding the characteristics of Zonotopes and p-Zonotopes can be found in the prior literature (Alanwar et al., 2019; Althoff et al., 2009).

For example, the 1D p-Zonotope shown in Figure 3(a) is formulated by considering the mean of distributions as lying within [−2, 2] and the covariance as lying between [0, 4]. Thereafter, we formulate the three parameters of a 1D p-Zonotope as follows: . Similarly, the 2D p-Zonotope shown in Figure 3(b) is formulated by considering the bounds on the mean and covariance of the first state to be [−2, 2] and [0, 4], respectively, and of the second state to be [−10, 10] and [0, 6], respectively. Mathematically, we represent this 2D p-Zonotope by and .

To perform set operations on stochastic reachable sets using p-Zonotopes, we require the following important properties (Althoff et al., 2009; Girard, 2005):

Minkowski sum of addition: The addition of two p-Zonotopes and is given by
4a
4b
where ⊕ is known as the Minkowski sum operator.
Linear map: Given a matrix and a p-Zonotope ,
4c
Translation: Given a vector and a p-Zonotope ,
4d
Order reduction: For any p-Zonotope with and a given maximum allowable order of e_max, the p-Zonotope after order reduction is given by
4e
where denotes the generator matrix after order reduction (Girard, 2005). To achieve this order reduction, the generators of G₁ are first sorted based on the Frobenius norm. Thereafter, one holds the first (e_max − 1) columns of G₁ which are denoted by G_>, and encloses the remaining columns of G₁, which are denoted by G_<, into a smallest box (interval hull), such that where rs(·) denotes the row-sum of a matrix.
γ -confidence: For with and Z = 〈c, G〉, transforming a p-Zonotope to a confidence zonotope, which is denoted by Z_γ, is given by
4f
where γ denotes the pre-defined parameter at the threshold of the tail of a p-Zonotope. This property is used to estimate the size of the p-Zonotope for a given γ. Note that, is well-known, where erf(·) denotes the error function (Oldham et al., 2008). Given l probabilistic generators in , the probability of lying outside the γ-confidence set is . When the state vector is considered as clock bias and clock drift, this variable r is analogous to the term timing risk associated with the γ-confidence set. As a simple example, for a predefined parameter γ = 3.5 and number of generators l = 2, the timing risk r = 1.9e⁻³. Extending this further, Section 4.4 of our proposed SR-DKF algorithm computes the timing risk associated with violating the pre-defined safety AL. For a more theoretical explanation regarding risk calculation for p-Zonotopes, we refer the readers to previous work (Althoff et al., 2009).

3 SYSTEM MODEL

Our system model comprises a network of power substations. Each power substation is equipped with a clock, PMU, GPS a receiver antenna, a communication transmit/receive antenna and a processing unit. In this work, for simplicity, we use the terms GPS receiver and power substation interchangeably. In the proposed SR-DKF algorithm, we mathematically represent the inter-connected network of power substations, each mounted with one GPS receiver, as an un-directed graph , where M indicates the set of nodes and E denotes the set of connections between the nodes. Each node represents a GPS receiver and each edge indicates a secure bidirectional communication link between the power substations. The total number of receivers in the network are |M|, where |·| denotes the cardinality of a set. Our SR-DKF algorithm does not require the network of power substations to be fully connected, hence,

The neighborhood set of any ith receiver ∀i ∈ M is represented by M_i and indicates the subset of nodes among the M nodes that can communicate with the ith receiver, including itself. Therefore, the associated cardinality of the neighborhood set is denoted by |M_i| where |M_i| ≥ 1 ∀i ∈ M. Every receiver in the network has prior knowledge regarding the pre-computed static position of all the receivers in its neighborhood set.

We define the true global clock time (which is expressed in GPS time) at the ith receiver and at the kth time iteration to be . Furthermore, based on prior work (Coleman & Beard, 2020; Khalife & Kassas, 2017), we define a local state vector at each ith receiver as , where cδt_i,k denotes the receiver clock bias, denotes the receiver clock drift and c represents the speed of light. The global clock time and the local state vector x_i,k are related as

where denotes the local measured time output (known) at the ith receiver. Note that, while the local measured time can be determined from the power substation at each time iteration, the global clock time is unknown and needs to be estimated. In practice, one way to estimate local measured time is to propagate the previously authenticated global time (e.g., when estimated timing risk from our proposed SR-DKF algorithm is less than the user-acceptable value) at a past time iteration over shorter time spans using independent clock model at each power substation or via communication links between base stations.

As explained in Section 1, the estimate of global clock time at each ith receiver is used by the PMU for time-tagging the recorded phasor values. Therefore, the accuracy of the global clock time depends on the accuracy of the estimated receiver’s clock bias cδt_i,k, which in-turn depends on the estimated clock drift . When a receiver is spoofed, the estimate of local state vector x_i,k is manipulated, thereby providing a misleading estimate of global clock time. Given this direct dependency, we investigate a distributed state estimation framework that considers x_i,k to be the state vector at the ith receiver.

The state-space model representing the true system at the ith receiver ∀i ∈ M is given by

with

where

– z_i,k denotes the 2N_i × 1 GPS measurement vector and N_i denotes the number of GPS satellites at the ith receiver. Here, and denote the received pseudorange and Doppler value from the nth satellite ∀n ∈ {1, ⋯, N_i}, respectively. Also, p_i and denote the known position and zero velocity of the ith static receiver, respectively, and pⁿ and denote the position and velocity of the nth satellite computed from the broadcast navigation message, respectively.
– F_i denotes the 2 × 2 state transition matrix, such that , where δt denotes the update rate.
– H_i denotes the 2N_i × 1 GPS measurement model, such that , where I_a denotes an a × a identity matrix and 0_a denotes an a × a zero matrix;
– ω_{i k} denotes the 2 × 1 process noise vector that includes all the clock-related errors that depict its short-term and long-term stability characteristics;
– ν_i,k denotes the 2N_i × 1 measurement noise vector that is GPS pseudoranges and Doppler values, caused by errors in satellite clock corrections and ionospheric and tropospheric effects, among others.

At any ith receiver, GPS measurements are received from all receivers in its neighborhood set M_i along with their estimated attack status, which quantifies their confidence in the authenticity of their received measurements. A detailed explanation regarding estimations of attack status is provided later in Section 4.3. In addition, each jth receiver j ∈ M_i sends the known local measured time to the ith receiver, which represents the time at which the message is transmitted.

Based on this, the state-space model representing the relationship between the global clock time of the ith receiver ∀i ∈ M and any jth receiver belonging to its neighborhood set ∀j ∈ M_i is given by

where D_ij,k denotes the average transmission delay in sending the data packet from any jth receiver at its local measured time to the ith receiver. The average transmission delay D_ij,k can be estimated from the two-way exchange of messages over time. Also, ω_ij,k denotes the zero-mean random term associated with the transmission delay between the ith and the jth receivers. Without loss of generality, we can consider the data packet from each jth receiver in the neighborhood set j ∈ M_i to be received at the ith receiver at its own local measured time with k_j ∈[k, k + 1). For detailed explanation regarding the clock modeling for distributed framework, the reader is referred to previous work (Bhamidipati et al., 2022; Mo et al., 2013)

From Equations (5) and (7), we formulate a true system that depicts the relationship between the clock biases at the ith receiver and the jth receiver ∀j ∈ M_i as

At any ith receiver, we formulate the relationship between the state vector x_i,k and the measurement model of the data sent from any jth receiver ∀j ∈ M_i as follows:

where denotes the measurement vector, such that , which can be observed from Equations (6), (7), and (8). Also, H_ij denotes the associated measurement model such that with N_j as the number of visible GPS satellites at the jth receiver and v_ij,k denotes the noise in the measurement vector z_ij,k that includes noise in the jth receiver GPS pseudorange and the transmission delay ω_ij,k. Therefore, the ith receiver constructs the measurement vector z_ij,k using the data packet sent from the jth receiver comprising and .

Intuitively, the measurement vector z_ij,k provides a noisy estimate of the ith receiver’s clock bias given the GPS measurements of the ith receiver. On the other hand, the measurement vector z_ij,k provides a noisy estimate of the ith receiver’s clock bias given the GPS measurements of the jth receiver. At any ith receiver, the measurement model H_ij depicts that only the GPS pseudoranges from the jth receiver are utilized in our proposed DKF formulation, while Doppler values are discarded. This is because the clock drift at the ith receiver is independent of the clock drift at any jth receiver; thus, we cannot utilize only one data packet sent from any jth receiver to the ith receiver to construct a noisy estimate of the ith receiver’s clock drift given the jth receiver’s GPS measurements. However, these additional measurements can be constructed by sending two (rather than one) consecutive data packets from the jth receiver to the ith receiver and thereafter by estimating the variation in the difference between the associated local measured times at both the receivers across the data packets. This procedure is explained in Section 3.2 of Mo et al., 2013 and is left as a potential future extension of our work.

For the system model described in Equations (5)–(9), we refer the readers to Appendix A for the mathematical formulation of the point-valued DKF comprising of two steps, i.e., time update and measurement update. The two point-valued DKF techniques (conventional and adaptive) are used as benchmarks for validating the proposed SR-DKF algorithm, which is discussed via simulated experiments in Section 5.

In this work, we feature several specific design considerations/assumptions. First, we consider the power substations as independent units that are installed in open-sky areas and geographically separated by large distances. Also, we make no assumption regarding whether the power substations in the network have the same hardware specifications, i.e., they can have GPS receivers of varying quality. Therefore, we do not consider the receiver clock drift across the network to be the same. We also do not account for correlations in measurement errors across receivers, i.e., we model them as pairwise independent. Given the open-sky conditions, we also assume that the multipath effects have a negligible impact on received measurements. Furthermore, our SR-DKF algorithm does not account for the presence of nominal orbit errors or Signal-In-Space (SIS) anomalies in the broadcast navigation message that can potentially affect multiple receivers and cause correlated errors. These assumptions are justified because there exists rich literature (Jiang & Groves, 2012; Ochieng et al., 2003; Xu et al., 2020) that addresses other sources of measurement error, including SIS anomalies, atmospheric effects, multipath, and others that can be independently added as a pre-processing step to the proposed SR-DKF algorithm.

4 THE PROPOSED SR-DKF: A SET-VALUED SECURE STATE ESTIMATION ALGORITHM

Based on the design aspects discussed in Section 1.1, we begin by outlining the architecture of the proposed SR-DKF algorithm and later describe the detials in a step-by-step fashion. While our current mathematical formulation considers only GPS measurements, the proposed SR-DKF can be easily generalized to incorporate multiple Global Navigation Satellite System (GNSS) constellations. This will require the system model to be updated to incorporate clock bias associated with all GNSS constellations to be considered with the state transition matrix and the measurement model modified accordingly. The proposed SR-DKF algorithm can also be updated to incorporate higher-order clock error models (Coleman & Beard, 2020) by considering additional states and modifying the state transition matrix accordingly.

Figure 4 presents a flowchart of the proposed SR-DKF algorithm that is executed independently at the ith receiver/power substation, and comprises four main modules, namely time update of set-valued DKF, spoofing attack mitigation, measurement update of set-valued DKF, and timing risk analysis.

FIGURE 4

Architecture of the proposed SR-DKF algorithm, which is executed independently at the ith receiver, comprises four main modules, namely time update of set-valued DKF, spoofing attack mitigation, measurement update of set-valued DKF, and timing risk analysis.

First, we perform a set-valued time update to compute the predicted p-Zonotope of the local state vector. Next, at each ith receiver, we independently analyze the received measurements from visible satellites in the set-valued domain by comparing them against the predicted p-Zonotope and known measurement error bounds in authentic conditions. Based on this, we estimate a metric termed as attack status that lies between 0 − 1, where 0 indicates that the received measurements comply with the known measurement error bounds in authentic conditions and thus are likely to be authentic whereas 1 indicates that the measurements are likely spoofed. In an intuitive sense, the attack status adaptively scales the known p-Zonotope of measurement error bounds in authentic conditions to account for the received measurements. Thereafter, the GPS measurements and their estimated attack status are asynchronously broadcast across receivers within the network. In the measurement update at ith receiver, the GPS measurements and their attack status from neighbors, i.e., j ∈ M_i, are collectively processed to compute the corrected p-Zonotope of the local state vector. The center mean of the corrected p-Zonotope along with the measured time at the ith receiver is used to estimate the point-valued GPS time that is later provided to the PMUs. The center uncertainty and covariance of the corrected p-Zonotope are further analyzed for the probability of their intersection with an unsafe set to compute the associated timing risk. In this context, we define unsafe set by a set of states, i.e., receiver clock bias and clock drift that violate the pre-defined safety AL that was defined earlier in Section 1. Thereafter, the process then repeats for next iteration.

In our proposed SR-DKF algorithm, we chose to perform spoofing attack mitigation; however, another potential approach includes attack detection/exclusion, provided there are sufficient number of non-spoofed measurements in the receiver network.

4.1 During Initialization: Pre-Defined Process and Measurement Error Bounds in Non-Spoofed Conditions

As explained earlier in Section 1 and from prior literature (Bhamidipati & Gao, 2018), it is justified to consider the error statistics to be known and bounded during non-spoofed (i.e., authentic) conditions. We can perform offline empirical analysis to calculate the process and measurement noise bounds for each GPS receiver under authentic conditions by analyzing the associated historical GPS data. Specifically, using the three-parameter representation of p-Zonotopes we can efficiently represent the different non-zero mean but bounded stochastic errors in a process and measurement model. Some examples of error sources were discussed earlier in Equation (6). Additional details regarding process and measurement error bounds considered in our experimental validation are provided in Section 5. For all receivers in the network, i.e., ∀i ∈ M, we model the process and measurement noise bounds as p-Zonotopes using Equation (1), which are represented by and , respectively, as follows:

where and denote the generator and covariance matrices of p-Zonotope that represents process noise bounds. Similarly, and denote the generator and covariance matrices of p-Zonotope that represents measurement noise bounds. If one considers the center of p-Zonotope for the process and measurement noise bounds to be a zero vector, the biases in errors associated with the empirical analysis are modeled via and , respectively. Furthermore, the covariance matrices and represent the characteristics of p-Zonotopes and are different from the point-valued process noise covariance, which is denoted by Q_i,k, and point-valued measurement covariance, which is denoted by R_i,k. More details regarding the matrices Q_i,k and R_i,k are provided in Appendix A.

4.2 Set-Valued DKF

We formulate the set-valued DKF by applying stochastic reachability to the point-valued DKF as described in Appendix A. Similar to point-valued DKF, the set-valued DKF also performs time and measurement updates at each ith receiver to compute the predicted and corrected stochastic reachable sets of the state vector x_i,k, respectively. To perform spoofing attack mitigation and timing risk analysis across a network of GPS receivers, we analyze the predicted and corrected stochastic reachable sets of the state estimation error, respectively. To represent the predicted and corrected point-valued state estimates at the kth iteration by and , respectively, the associated estimation error is given by and , where and .

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 and , respectively, are given by

where and denote the generator and covariance matrices of the corrected p-Zonotope of the state estimation error. The corrected p-Zonotope is estimated during the measurement update of set-valued DKF, which is explained later in Section 4.2.2. Similarly, and 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 4.2.1. Note that the aforementioned covariance matrices and represent the Gaussian characteristics of the p-Zonotopes and are different from the predicted point-valued state covariance matrix, which is denoted by , and the corrected point-valued state covariance matrix, which is denoted by . More details regarding the matrices and are provided in Appendix A.

At the kth iteration, we utilize the translation set property from Equation (4d) to represent the predicted and corrected p-Zonotopes of the local state vector and , respectively, as and . Therefore, the predicted and corrected p-Zonotopes of the local state vector are given by Equation (12a) and Equation (12b), respectively. Note that translating the center of p-Zonotope as seen in Equation (12) 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_i,k and state estimation error Δx_i,k are the same. From Equations (12a)–(12b) and prior literature (Althoff et al., 2009; Shi et al., 2015), it has been established that the center mean of predicted and corrected p-Zonotopes is equivalent to their point-valued DKF estimates.

12a

12b

4.2.1 Predicted p-Zonotope of the Local State Vector via Time Update

At the kth iteration, we perform time update of set-valued DKF independently at each ith receiver in the network ∀i ∈ M to compute the predicted p-Zonotope of the local state vector. As explained earlier in Equation (12a), the center mean of the predicted p-Zonotope of the local state vector is the same as the point-valued estimate of the conventional DKF, whose mathematical formulation is described in Appendix A. Thereafter, we estimate the generator and covariance matrices of the predicted p-Zonotope of the state estimation error by utilizing the corrected p-Zonotope of the previous (k – 1) th iteration and the p-Zonotope of the process noise bound. To execute these steps, we first formulate the point-valued state estimation error as

Given that are independent random variables, we then apply the properties of p-Zonotopes discussed earlier in Equations (4a)–(4b) for converting the point-valued representation of Equation (13) to estimate the predicted p-Zonotope as seen in Equations (14a)–(14b). We refer readers to relevant literature (Althoff & Dolan, 2014; Althoff et al., 2009) for examples of the conversion of point-valued state space Equations to set-valued formulations.

14a

such that

14b

For the example shown in Figure 5, the predicted p-Zonotope obtained after the Minkowski sum in Equation (14a) is indicated by the Parula colormap. The predicted p-Zonotope of the state estimation error is used to perform spoofing attack mitigation as described in Section 4.3.

FIGURE 5

Visualization of the predicted p-Zonotope of the state estimation error after the time update of set-valued DKF. The 2D p-Zonotope is sliced to be shown as two 1D p-Zonotopes, one each for (a) Set-valued time update: clock bias and (b) Set-valued time update: clock drift. The Minkowski sum of the pre-defined p-Zonotope of the process noise, which is indicated in gray, and the linear map of corrected p-Zonotope at the (k – 1)th iteration, indicated in orange, provide an estimate of the predicted p-Zonotope for the kth iteration. The predicted p-Zonotope is indicated by the Parula colormap (Nuñez et al., 2018), where the range of colors from blue to yellow denotes the lowest to highest values of probability represented by this p-Zonotope.

4.2.2 Corrected p-Zonotope of the State Estimation Error via Measurement Update

In the measurement update executed at ith receiver, we compute the corrected p-Zonotope of the local state vector. As explained earlier in Equation (12b), the center mean of the corrected p-Zonotope of the local state vector is the same as the point-valued estimate of the conventional DKF, whose mathematical formulation is described in Appendix A. Thereafter, we estimate the generator and covariance matrices of the corrected p-Zonotope of the state estimation error by utilizing the p-Zonotopes of measurement noise bounds associated with all the receivers in the neighborhood set M_i and predicted p-Zonotope of the state estimation error. In our distributed framework, as explained earlier in Section 3, each ith receiver utilizes both pseudoranges and Doppler values received at its location but utilizes only the pseudoranges received at other receivers in the neighborhood set (M_i – i).

Given that we only know the measurement error bounds under authentic conditions, we designed an adaptive framework to perform the measurement update of the set-valued DKF. The adaptive gain matrix associated with ith receiver measurement is denoted by and that associated with received measurements from neighborhood set j ∈ (M_i – i) is denoted by . Further details regarding the adaptive gain matrix formulation are provided in the spoofing attack mitigation module, which is explained later in Section 4.3. To execute this function, we first formulate the point-valued state estimation error as

where z_i,k and H_i denotes the received measurements from ith receiver in the network and the associated measurement model, respectively, as defined earlier in Equation (6). As described in Appendix A, z_ij,k and H_ij,k denote the received measurements from other receivers in the neighborhood set of the ith receiver and the associated measurement model, respectively, as defined earlier in Equation (9).

Similar to the time update of set-valued DKF explained earlier in Section 4.2.1, we observe that are independent random variables. Based on this observation, we the apply the properties of p-Zonotopes discussed earlier in Equations (4a)–(4b) for converting the point-valued representation of Equation (15) to set-valued stochastic reachable sets as seen in Equation (16a).

16a

We further simplify using Equation (4c) to obtain Equation (16b). Intuitively, the estimated generator and covariance matrices, i.e., and , provide bounds on the unknown state estimation error associated with the point-valued state estimate .

16b

In the example shown in Figure 6, the corrected p-Zonotope obtained after the Minkowski sum in Equation (16a) is indicated by a Parula colormap. Also, given that the size of the generator matrix shown in Equation (16b) is increasing during each measurement update, we apply the set property of p-Zonotope discussed earlier in Equation (4e) to reduce the order of the generator matrix of the corrected p-Zonotope of the state estimation error . The corrected p-Zonotope of the state estimation error is analyzed to perform timing risk analysis as described in Section 4.4.

FIGURE 6

Visualization of the corrected p-Zonotope of the state estimation error after the measurement update of set-valued DKF. The 2D p-Zonotope is sliced to be shown as two 1D p-Zonotopes, one each for (a) Set-valued measurement update: clock bias and (b) Set-valued measurement update: clock drift. The Minkowski sum of the linearly-mapped predicted p-Zonotope at the kth iteration, indicated in red, and the pre-computed measurement noise bounds, indicated in gray, provide an estimate of the corrected p-Zonotope for the kth iteration. The corrected p-Zonotope is indicated by the Parula colormap (Nuñez et al., 2018), where the range of colors from blue to yellow denotes the lowest to highest values of probability represented by this p-Zonotope.

4.3 Spoofing attack mitigation and data transfer

Before performing measurement updates at the ith receiver and kth iteration, we independently utilized the predicted p-Zonotope of the state estimation error from Equation (14) and the pre-computed measurement noise bound to estimate the p-Zonotope of expected measurement innovation. For this application, we utilized the point-valued measurement innovation given by . More details are discussed in Appendix A. We then apply the set properties of p-Zonotopes shown in Equations (4a) and (4c) on the point-valued measurement innovation ϵ_i,k to compute the p-Zonotope of expected measurement innovation, which is denoted by , as follows:

The p-Zonotope of expected measurement innovation maps each point-valued measurement innovation ϵ_i,k to a scalar value that is denoted by α_i,k. Based on this, at each ith receiver, we independently compute a metric termed as attack status by normalizing the obtained value α_ij,k with the value corresponding to the center mean of the p-Zonotope and subtracting this value from one. Thus, the estimated attack status lies between [0, 1]. As seen earlier in Equation (10), the pre-computed process and measurement noise bounds are defined for non-spoofed conditions. Therefore, in a non-spoofed condition, such as that shown in Figure 7(a), the point-valued measurement innovation ϵ_i,k and its p-Zonotope of expected measurement innovation are in close agreement, and low attack status is obtained. However, during spoofing, as shown in Figure 7(b), the point-valued measurement innovation does not comply with this estimated p-Zonotope, and thus a high value of is observed. Note that our aforementioned analysis on spoofing attack mitigation is conducted in vector domain. Therefore, even when a subset of visible GPS satellites is affected by spoofing, the estimated attack status remains high at , indicating low overall confidence in the ith receiver.

FIGURE 7

Visualization of the attack status during (a) non-spoofed conditions when attack status ; (b) spoofed conditions when attack status . The known p-Zonotope of measurement error bound before and after scaling are represented by a Parula colormap (Nuñez et al., 2018) where the range of colors from blue to yellow denotes the lowest to highest values of probability represented by this p-Zonotope. By adaptively scaling the known p-Zonotope of measurement error bounds in authentic conditions based on the estimated attack status, the received measurement innovation is efficiently captured by the center uncertainty (i.e., the generator matrix).

Each ith receiver in the network i ∈ M broadcasts its received GPS measurement z_ij,k as was defined earlier in Equation (15), as well as the estimated attack status to all the receivers in the neighborhood set j ∈ (M_i – i). The set (M_i – i) simply implies that the ith receiver need not broadcast GPS measurements to itself. Similarly, each ith receiver also receives the GPS measurements z_ij,k and their associated attack status from one or more receivers in its neighborhood set j ∈ (M_i – i). In an intuitive sense, the attack status obtained from each jth receiver in the neighborhood set indicates its belief/confidence in the received GPS measurement z_ij,k. Thereafter, at the ith receiver, the attack status is used to compute the adaptive gain matrices and as follows:

where R_ij,k denotes the measurement covariance associated with the received measurements z_ij,k at the ith receiver from jth receiver ∀j ∈ (M_i – i). The mathematical formulation of z_ij,k, R_ij,k and H_ij,k were described earlier in Equation (15). More details regarding the mathematical formulation of the point-valued state covariance matrix can be found in Appendix A. Formulated adaptive gain matrices and are utilized to perform the set-valued measurement update of the SR-DKF algorithm, as was described earlier in Equation (16).

By comparing Equation (18) and Equation (16), the adaptive gain matrices scale the known p-Zonotope of measurement error bounds in authentic conditions adjusts to comply with the received measurements. This is reflected in Figure 7(a) that represents non-spoofed conditions, wherein the adaptive scaled p-Zonotope remains close to its original size (i.e., before scaling) with the received measurement innovation captured by the center uncertainty (i.e., generator matrix). Similarly, in Figure 7(b) which represents the spoofed conditions, the size of the adaptively-scaled p-Zonotope increases significantly with the center uncertainty accounting for the corresponding received measurement innovation. Note that, unlike the adaptive point-valued DKF techniques discussed earlier in Section 1 that accounts for stochastic biases as a point-valued covariance, the proposed SR-DKF algorithm models stochastic biases in set-valued domains using center uncertainty (via the generator matrix).

4.4 Timing Risk Analysis

We evaluate the probability that the corrected p-Zonotope of state estimation error derived in Equation (16) falls within an unsafe set, i.e., a set of states that violate the pre-defined safety AL = ±1 μs. As shown in both Figures 8(a) and 8(b), the unsafe set, which is denoted by , is represented by two gray strips, i.e., and .

FIGURE 8

Visualization of the corrected p-Zonotope of state estimation error and its intersection with the unsafe set . (a) Over-approximation of corrected p-Zonotope: shows an example that over-approximates the corrected p-Zonotope by four polytopes, which are denoted by D₁ to D₄; (b) Top view of finite polytopes: shows the top-view of the finite polytopes and an example of the intersection region of D₃ polytope with unsafe set , which is indicated by red cross marks.

We perform timing risk analysis by over-approximating the corrected p-Zonotope using a finite series of polytopes (Althoff et al., 2009). To over-approximate the corrected p-Zonotope efficiently, we first apply the γ-confidence set property of p-Zonotope that was discussed earlier in Equation (4f) to compute a series of (L +1) confidence zonotopes, which are denoted by , for the corresponding pre-defined values of γ = m₀, m₁, m₂, ⋯, m_L where m₀ > m₁ > m₂ > ⋯ > m_L. Note that the confidence zonotope Z_m₀ represents the upper bound at the threshold of the tail of the corrected p-Zonotope , and therefore, the larger the value of m₀, the more exact the representation of the corrected p-Zonotope by a confidence zonotope will be. Furthermore, as the value of L increases, the finite series of polytopes more closely approximate the corrected p-Zonotope, but at the expense of an increased computational cost.

Using the (L +1) confidence zonotopes Z_m₀ to Z_{m_L}, we formulate a series of L polytopes, which are denoted by D₁ to D_L, as follows: , where the set operation for any given sets S_a and S_b. Further details regarding the over-approximation of p-Zonotopes by a series of polytopes are discussed in Althoff et al., 2009. An example of the over-approximation of p-Zonotope by four polytopes D₂ to D₄ is shown in Figure 8(a). We then compute the area of intersection of each lth polytope D_l ∀l ∈ {1, ⋯, L} with the unsafe set . Figure 8(b) shows the top-view of an over-approximated series of polytopes and an example of the intersection region of D₃ polytope with the unsafe set .

In the context of p-Zonotopes (Althoff et al., 2009; Georgescu et al., 2019), the timing risk at ith receiver, which is denoted by r_i,k, is theoretically defined as

For computational tractability, we re-wrote the timing risk r_i,k, which is seen in Equation (20), as a summation of two terms. The first term represents the risk associated with the confidence zonotope Z_m₀ as exceeding the unsafe set and is formulated based on the γ-confidence property of p-Zonotope described earlier in Equation (4f). The second term denotes the sum of the intersection area across polytopes D_l, which is denoted by , with the unsafe set times the maximum probability of D_l th polytope, which is denoted by p_max,l. The maximum probability of D_l th polytope is at its top face and equals the probability of the corrected p-Zonotope at that level. Note that, since the area under the p-Zonotope does not equal one, the associated timing risk r_i,k described in Equation (20) also does not necessarily need to be less than one. Therefore, based on the theoretical formulation explained in Althoff et al., 2009, if r_i,k < 1, the estimated timing exceeds the safety AL with a risk of r_i,k whereas if r_i,k ≥ 1 then the estimated timing exceeds the safety AL with a risk of 1.

4.5 Implementation Details

Algorithm 1 summarizes the implementation steps of the proposed SR-DKF algorithm described earlier in Sections 4.2, 4.3, and 4.4. To compute the adaptive gain matrix at each kth iteration, we also propagate the point-valued state covariance matrix with the predicted estimate denoted by and the corrected estimate denoted by . More details regarding the formulation of predicted and corrected point-valued state covariance can be found in Appendix A. During initialization, i.e., at iteration k = −1, we denote the initial point-valued mean and covariance of the state vector by and , respectively. Furthermore, we consider all receivers to be authentic during initialization. Based on this, we represent the initial p-Zonotope of the state estimation error to be zero-mean center, with a non-zero generator and over-bounding covariance matrices. For simplicity and without the loss of generality, we process the measurements from different receivers in a single neighborhood set simultaneously. Alternatively, for a more generalized framework, one can process the measurements in an asynchronous manner, i.e., execute the proposed SR-DKF algorithm at every time sub-iteration when measurements are received from any jth receiver in the neighborhood set.

ALGORITHM 1

Proposed SR-DKF algorithm for spoofing attack mitigation and timing risk analysis in networked power systems

5 RESULTS AND ANALYSIS

We validate the robustness of the proposed SR-DKF algorithm to not only mitigate the effects of spoofing but also to estimate the secure GPS timing and its associated timing risk across all the receivers in a network. We perform set operations of stochastic reachability and transform across various set representations, including polytopes, zonotopes, p-Zonotopes, and others, using a publicly-available MATLAB toolbox, known as the COntinuous Reachability Analyzer (CORA) (Althoff, 2016).

As shown in Figure 9, we designed a simulated network of seven GPS receivers, |M| = 7, that are spatially distributed across the 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-defined GPS simulator known as GPS-SIM-SDR (Bhamidipati & Gao, 2019; Ebinuma, 2018) to simulate the raw GPS signals received at each location in the network. We then post-processed the simulated GPS signals using a MATLAB-based software-defined radio known as SoftGNSS (Paul, 2011). In our simulation setup, we set the average transmission delay to a pre-surveyed/constant value, which is also considered in modeling the measurement vector z_ij,k in the system model of our proposed SR-DKF algorithm.

FIGURE 9

Geographical map showcasing a simulated network of seven spatially distributed GPS receivers across the US

We defined the simulated errors observed in the system during authentic conditions as having the following characteristics:

For the process noise in local state vector (the first state cδt_i,k), the mean and covariance of the time-varying Gaussian distribution lies between [−36 m, 36 m] and [0 m, 60 m], respectively; similarly, for the process noise in its clock drift (the second state ), the time-varying mean and covariance lies between [−0.024 m/s, 0.024 m/s] and [0 m/s, 0.03 m/s], respectively.
For the measurement noise in GPS pseudoranges, the mean and covariance of the time-varying Gaussian distribution lies between [−2.5 m, 2.5 m] and [0 m, 3.68 m], respectively; similarly, for the measurement noise in GPS Doppler, the time-varying mean and covariance lies between [−0.33 Hz, 0.33 Hz] and [0 Hz, 0.05 Hz], respectively. We consider all the receivers in this simulated setup to have the same error bounds. We refer to prior studies (Hetet, 2000; Kuusniemi et al., 2004) for information on deriving the empirical measurement noise covariance of pseudoranges and Doppler values as a function of carrier-to-noise density (C/N0). These values are formulated by setting C / N 0 = 38 dB-Hz to depict authentic conditions.
For the initial state error bounds of local state vector (the first state cδt_i,k), the mean and covariance of the time-varying Gaussian distribution lies between [−0.01 μs, 0.01 μs] and [0 μs, 0.04 μs], respectively; similarly, for the initial state error bounds in the second state , the time-varying mean and covariance lies between [−0.02 ns/s, 0.02 ns/s] and [0 ns/s, 0.06 ns/s], respectively.

As discussed earlier in Algorithm 1, we initialized our SR-DKF algorithm by pre-defining the following p-Zonotopes ∀i ∈ M including process noise , measurement noise , and initial state estimation error . We formulate these p-Zonotopes from the aforementioned errors bounds in the same manner as the 2D example discussed in Section 2. Additional complexities associated with the simulation setup and exhaustive sensitivity analysis for different ranges of hyperparameters will be explored in our future work.

5.1 During a Coordinated Signal-Level Spoofing Attack

In the first set of experiments, we consider sparse connectivity among the seven receivers in the network. Figure 10(a) showcases the connectivity matrix where a green tick indicates the presence of a bidirectional communication link while a red cross indicates otherwise. Note that each receiver is connected to itself by default.

FIGURE 10

Performance comparison of the three state estimation techniques. a) sparse connectivity across the network of receivers; b) set-valued state estimation using the proposed SR-DKF algorithm; c) point-valued adaptive DKF; d) point-valued conventional DKF; In (b)-(d) the black-dotted vertical lines indicate the start and stop times of spoofing attacks while the magenta-dotted horizontal lines indicate the threshold of ±26.5 μs set by the IEEE C37.118.1a-2014 standards, as a reference.

We induce a simulated coordinated signal-level spoofing attack in the simulated GPS measurements received at two GPS receivers in the network, i.e., Rx : 1, which is located in Stanford, CA and Rx : 5, which is located in Austin, TX. Between the time duration k = 40 −1040 s, the first spoofing attack manipulates the GPS time at Rx : 5 to deviate at a rate of 100 ns/s. The second spoofing attack, which occurs between k = 800 −1300 s, deviates the GPS time of Rx : 1 at a rate of 400 ns/s.

We validated the performance of the proposed SR-DKF algorithm against two benchmark point-valued state estimation techniques, namely adaptive DKF, which is shown in Figure 10(c), and conventional DKF, which is shown in Figure 10(d). The mathematical formulations for these point-valued techniques are provided in Appendix A; the measurement noise covariance R_i,k is estimated by setting the forgetting factor at ψ = 0.3 in the case of adaptive DKF and ψ = 1 in the case of conventional DKF. The quantitative maximum error values of the receiver network for the aforementioned state estimation approaches are listed in Table 1. The conventional DKF shown in Figure 10(d) exhibits high state estimation errors in the receivers that are being spoofed, i.e., Rx : 5 and Rx : 1, with maximum clock bias errors of 69.24 μs and 13.77 μs, respectively. The conventional DKF also corrupts the timing at other non-spoofed receivers in the network, e.g., Rx : 3 and Rx : 2 exhibit large state estimation errors in clock bias, whose maximum values are 100.07 μs and 42.88 μs, respectively. Furthermore, while the adaptive DKF, which is shown in Figure 10(c), demonstrates higher state estimation accuracy than the conventional DKF, this point-valued technique still affects the receivers that are not directly spoofed, e.g., Rx : 3 exhibits state estimation errors with a maximum value of 2.47 μs. Based on the maximum error values listed in Table 1, we observe that both point-valued state estimation techniques violate the pre-defined safety AL = ±1 μs. By contrast, our proposed SR-DKF algorithm, which is shown in Figure 10(b), provides a secure estimate of GPS timing with a maximum clock bias error of 0.42 μs and clock drift error of 2.23 ns/s across all receivers; thus, the SR-DKF algorithm complies with the acceptable timing standards that are commonly practiced in the power community.

View this table:

TABLE 1

Comparison among the maximum error values of the state vector estimated via three approaches: the proposed SR-DKF algorithm, point-valued adaptive DKF, and point-valued conventional DKF

In Figure 11, we plot the attack status of the GPS measurements z_i,k at each receiver in the network. In accordance with Section 4.3, we demonstrate that the proposed SR-DKF algorithm successfully detects and mitigates spoofing in the GPS measurements reported by Rx : 1 and 5 while accurately categorizing other receivers, i.e., Rx : 2, 3, 4, 6, 7 as non-spoofed with attack status at all times. We also observe that, irrespective of the drift rate of spoofing, i.e., 400 ns/s in Rx : 1 and 100 ns/s in Rx :5, both attacks were quickly mitigated, thereby ensuring high robustness. Even after the spoofing ends, i.e., at k = 1040 s in Rx :5 and at k = 1300 s in Rx : 1, the GPS measurements still exhibit errors as their tracking loops have not yet converged to the authentic satellite signals. This effect is also captured by our proposed SR-DKF algorithm which has successfully assigned an attack status of 1 for both Rx :5 and Rx : 1.

FIGURE 11

Attack status of each receiver in the network as estimated by the proposed SR-DKF algorithm. This is computed by evaluating the point-valued measurement innovation against the estimated p-Zonotope of expected measurement innovation, as seen in Equation (17). Our SR-DKF quickly mitigates spoofing in Rx : 1, 5 and successfully estimates that other receivers have not been spoofed.

In Figure 12, we showcase the timing risk associated with the estimated GPS timing, i.e., the probability of estimated clock bias that violates the pre-defined safety AL = ±1 μs. This is computed by analyzing the intersection of the corrected p-Zonotope with the unsafe set , as discussed in Section 4.4. An increase or jump in measure of timing risk is based on the adaptive gain matrix that in turn depends on the estimated attack status of the receiver. By comparing results shown in Figure 10(a) and Figure 12, we observe a direct correlation between the value of timing risk to the number of reliable measurements available for the calculation of the 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., ≈ 10⁻⁷, for the entire duration of the experiment. By contrast, Rx : 3 has a highly fluctuating timing risk because two of its three neighbors have been spoofed. Similarly, when only one of three neighbors at Rx : 1 is spoofed, the order of timing risk magnitude is 10⁻³; by contrast, when only one of two neighbors are spoofed at Rx :5, the timing risk exhibits a higher order magnitude at 0.25. Therefore, our SR-DKF algorithm estimates a robust measure of the timing risk across all receivers in the network.

FIGURE 12

Robust timing risk analysis for each receiver in the network via the proposed SR-DKF algorithm. This is computed by evaluating the intersection of the corrected p-Zonotope of the state estimation error with the unsafe set as shown in Equation (20). Since risk requirements are not currently set by the power community, we pre-defined the safety AL = ±1 μs. Our proposed SR-DKF algorithm successfully quantifies the risk of violating the pre-defined safety AL. Note that the y-axis legend are not in the same range across all the receivers in the network. This was done to show the zoomed variation in the timing risk at each receiver under coordinated spoofing.

5.2 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 receivers in the neighborhood set |M| varies. Considering a total experiment duration of 100 s and the network of seven receivers shown in Figure 9, we induced a simulated meaconing attack between k = 10 −100 s in the GPS measurements that were received at Rx : 1 located at Stanford, CA.

In Figure 13, we analyzed the variation in the timing risk associated with the estimated GPS time of Rx : 1 as the number of receivers in the network M increased. For each receiver network size that ranges from |M| = 2 to |M| = 7, we consider the presence of a fully-connected bidirectional communication network and perform 50 Monte-Carlo runs for four meaconing magnitudes of 30 μs, 45 μs, 60 μs, and 100 μs. We successfully demonstrated that, as the number of receivers interacting with the Rx :1, i.e., |M₁|, increases, the associated timing risk decreases. Particularly for |M₁| > 5, the decrease in timing risk is on the order of ≤ 10⁻⁵ and hence considered insignificant. Similarly, the effect of meaconing magnitude on the timing risk decreases as the number of interacting receivers associated with Rx : 1 increases. Therefore, the proposed SR-DKF algorithm 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-connected network. The colors red, blue, cyan, and magenta indicate a meaconing attack of 30 μs, 45 μs, 60 μs, and 100 μs, respectively. We observe that the variation in timing risk decreases with an increase in the number of interacting receivers with Rx :1 and the magnitude of meaconing attack.

6 CONCLUSIONS

In summary, we developed a Stochastic Reachability-Based Distributed Kalman Filter (SR-DKF) algorithm to perform set-valued state estimation using a network of GPS receivers that estimate not only the point-valued mean but also the stochastic reachable set of the local state vector (i.e., receiver clock bias and receiver clock drift). By parameterizing the stochastic reachable set via probabilistic Zonotopes (p-Zonotopes), we derived both time and measurement update steps of the set-valued DKF to estimate the predicted and corrected p-Zonotope of 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 stochastic set (represented by the p-Zonotope). Second, we performed timing risk analysis at state-level via the probably of an intersection of the corrected p-Zonotope with an unsafe set, i.e., a set that violates the pre-defined safety alert limit (AL). To the best of our knowledge, no prior literature exists that addresses secure state estimation via stochastic reachability nor on the use of timing risk quantification during spoofing.

Under a simulated coordinated signal-level spoofing that affects two members of a network of seven sparsely-connected GPS receivers, 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 0.42 μs, both the point-valued state estimation techniques, i.e., conventional DKF and adaptive DKF, showed high estimation errors in clock bias that violate the pre-defined safety AL = ±1 μs, which is the commonly accepted timing accuracy for PMUs. Furthermore, by varying the number of receivers and attack magnitudes, we performed Monte Carlo runs to validate our finding that only a handful of GPS receivers (approximately five) are necessary to estimate GPS time with a low timing risk of 10⁻⁵.

HOW TO CITE THIS ARTICLE

Bhamidipati, S., & Gao, G. (2023). GPS spoofing mitigation and timing risk analysis in networked phasor measurement units via stochastic reachability. NAVIGATION, 70(3). https://doi.org/10.33012/navi.574

ACKNOWLEDGEMENTS

We would like to thank Ashwin Kanhere, Siddharth Tanwar, and Tara Mina for their help in revising the paper. We would also like to thank all members of the NAV Lab at Stanford University for their thoughtful input on this work.

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 United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

APPENDIX A PRELIMINARY FINDINGS OF THE POINT-VALUED DKF

Based on Talebi and Werner, 2018, the general framework of a point-valued DKF implementation performed independently at the ith receiver is shown in Equations (A.1)–(A.3). With as the expectation operator, we denote the point-valued mean and covariance during initialization, i.e., k = −1, by and , respectively, which are defined as follows:

A.1

The process and measurement noise at each ith receiver are modeled via a zero-mean Gaussian distribution as and , respectively, where Q_i,k and R_i,k represent the process and measurement noise covariance matrix, respectively. After initialization, the conventional point-valued DKF recursively executes two steps, including the time update and the measurement update.

During the time update, as seen in Equation (A.2), the corrected variables and from the previous time, i.e., at the (k – 1)th iteration, are propagated forward-in-time via a state transition matrix F_i to predict the estimates for the next kth iteration, which are given by and , respectively.

A.2a

A.2b

At the kth iteration, each ith receiver receives the GPS measurements and measurement noise covariance matrix 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 measurements are weighted via an estimated Kalman gain to obtain the corrected estimates of the point-valued mean and covariance of the state vector, i.e., and , respectively. In this context, z_i,k and H_i denotes the received measurements from the ith receiver itself and the associated measurement model, respectively, while z_ij,k, H_ij,k and R_ij,k denotes the measurements received from other receivers in the neighborhood set of the ith receiver and the associated measurement model and measurement noise covariance matrix, respectively. The optimal Kalman gain denoted by K_i,k and K_ij,k ∀j ∈ (M_i – i) is shown in Equation (A.3).

A.3a

A.3b

where gain matrix and . The mathematical formulation of H_i, z_i,k and for this work is discussed in Equation (5) and Equation (6) of Section 3. Similarly, the mathematical formulation of H_ij, z_j,k and is discussed in Equations (7), (8) and (9) of Section 3.

For an adaptive implementation of point-valued DKF, additional parameters, namely, forgetting factor ψ and pre-measurement residuals ϵ_i,k and ϵ_ij,k, are considered to provide an adaptive estimate the measurement noise covariance matrix R_i,k and R_ij,k at each time iteration. Here, with and with . Of note, we obtain the conventional point-valued DKF by setting ψ= 1.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

↵
1. Alanwar, A.,
2. Said, H., &
3. Althoff, M.
(2019). Distributed secure state estimation using diffusion kalman filters and reachability analysis. 2019 IEEE 58th Conference on Decision and Control (CDC). https://doi.org/10.1109/cdc40024.2019.9029929
↵
1. Al-Hammouri, A. T.,
2. Nordstrom, L.,
3. Chenine, M.,
4. Vanfretti, L.,
5. Honeth, N., &
6. Leelaruji, R.
(2012). Virtualization of synchronized phasor measurement units within real-time simulators for smart grid applications. 2012 IEEE Power and Energy Society General Meeting. https://doi.org/10.1109/pesgm.2012.6344949
↵
1. Allnutt, J.,
2. Anand, D.,
3. Arnold, D.,
4. Goldstein, A.,
5. Li-Baboud, Y.-S.,
6. Martin, A.,
7. Nguyen, C.,
8. Noseworthy, R.,
9. Subramaniam, R., &
10. Weiss, M.
(2017). Timing challenges in the smart grid (tech. rep.). National Institute of Standards; Technology. https://doi.org/10.6028/nist.sp.1500-08
↵
1. Althoff, M.
(2016). CORA 2016 manual. TU Munich, 85748.
↵
1. Althoff, M., &
2. Dolan, J. M.
(2014). Online verification of automated road vehicles using reachability analysis. IEEE Transactions on Robotics, 30(4), 903–918. https://doi.org/10.1109/tro.2014.2312453
↵
1. Althoff, M.,
2. Stursberg, O., &
3. Buss, M.
(2009). Safety assessment for stochastic linear systems using enclosing hulls of probability density functions. 2009 European Control Conference (ECC). https://doi.org/10.23919/ecc.2009.7074473
↵
1. Bhamidipati, S., &
2. Gao, G. X.
(2018). Multiple GPS fault detection and isolation using a graph-SLAM framework. Proceedings of the 31st International Technical Meeting of The Satellite Division of the Institute of Navigation (ION GNSS+ 2018). https://doi.org/10.33012/2018.16030
↵
1. Bhamidipati, S., &
2. Gao, G. X.
(2019). GPS multireceiver joint direct time estimation and spoofer localization. IEEE Transactions on Aerospace and Electronic Systems, 55(4), 1907–1919. https://doi.org/10.1109/taes.2018.2879532
↵
1. Bhamidipati, S., &
2. Gao, G. X.
(2020). GPS spoofing mitigation and timing risk analysis in networked PMUs via stochastic reachability. Proceedings of the 33rd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2020). https://doi.org/10.33012/2020.17757
↵
1. Bhamidipati, S.,
2. Iiyama, K.,
3. Mina, T., &
4. Gao, G.
(2022). Time-transfer from terrestrial gps for distributed lunar surface communication networks. IEEE Aerospace Conference.
↵
1. Bhamidipati, S.,
2. Mina, T. Y., &
3. Gao, G. X.
(2018). GPS time authentication against spoofing via a network of receivers for power systems. 2018 IEEE/ION Position, Location and Navigation Symposium (PLANS). https://doi.org/10.1109/plans.2018.8373542
↵
1. Bujorianu, L. M.
(2012). Applications of Stochastic Reachability. Stochastic Reachability Analysis of Hybrid Systems, 203–207. https://doi.org/10.1007/978-1-4471-2795-6_11
↵
1. Coleman, M. J., &
2. Beard, R. L.
(2020). Autonomous clock ensemble algorithm for GNSS applications. NAVIGATION, 67(2), 333–346. https://doi.org/10.1002/navi.366
↵
1. Crassidis, J. L., &
2. Junkins, J. L.
(2011). Optimal estimation of dynamic systems. Chapman; Hall/CRC. https://doi.org/10.1201/b11154
↵
1. Dagle, J.
(2018). Importance of synchrophasor technology in managing the grid. In Power electronics and power systems (pp. 1–11). Springer International Publishing. https://doi.org/10.1007/978-3-319-89378-5_1
↵
1. Ebinuma, T.
(2018). GPS-SDR-SIM.
↵
1. Georgescu, A.,
2. Gheorghe, A. V.,
3. Piso, M.-I., &
4. Katina, P. F.
(2019). Critical space infrastructures. In Critical space infrastructures (pp. 21–36). Springer International Publishing. https://doi.org/10.1007/978-3-030-12604-9_2
↵
1. Girard, A.
(2005). Reachability of uncertain linear systems using zonotopes. In Hybrid systems: Computation and control (pp. 291–305). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-540-31954-2_19
↵
1. Günther, C.
(2014). A survey of spoofing and counter-measures. Navigation, 61(3), 159–177. https://doi.org/10.1002/navi.65
↵
1. Heng, L.,
2. Work, D. B., &
3. Gao, G. X.
(2015). GPS signal authentication from cooperative peers. IEEE Transactions on Intelligent Transportation Systems, 16(4), 1794–1805. https://doi.org/10.1109/tits.2014.2372000
↵
1. Hetet, S.
(2000). Signal-to-noise ratio effects on the quality of gps observations. Fredericton, New Brunswick: Department of Geodesy and Geomatics, University of New Brunswick.
↵
1. Hofmann-Wellenhof, B.,
2. Lichtenegger, H., &
3. Collins, J.
(1992). Global positioning system. Springer Vienna. https://doi.org/10.1007/978-3-7091-5126-6
↵
1. Jiang, Z., &
2. Groves, P. D.
(2012). NLOS GPS signal detection using a dual-polarisation antenna. GPS Solutions, 18(1), 15–26. https://doi.org/10.1007/s10291-012-0305-5
↵
1. Khalajmehrabadi, A.,
2. Gatsis, N.,
3. Akopian, D., &
4. Taha, A. F.
(2018). Real-time rejection and mitigation of time synchronization attacks on the global positioning system. IEEE Transactions on Industrial Electronics, 65(8), 6425–6435. https://doi.org/10.1109/tie.2017.2787581
↵
1. Khalife, J., &
2. Kassas, Z. M.
(2017). Modeling and analysis of sector clock bias mismatch for navigation with cellular signals. 2017 American Control Conference (ACC). https://doi.org/10.23919/acc.2017.7963500
↵
1. Kuusniemi, H.,
2. Lachapelle, G., &
3. Takala, J. H.
(2004). Position and velocity reliability testing in degraded GPS signal environments. GPS Solutions, 8(4), 226–237. https://doi.org/10.1007/s10291-004-0113-7
↵
1. Li, W.,
2. Gong, D.,
3. Liu, M.,
4. Chen, J., &
5. Duan, D.
(2013). Adaptive robust kalman filter for relative navigation using global position system. IET Radar, Sonar & Navigation, 7(5), 471–479. https://doi.org/10.1049/iet-rsn.2012.0170
↵
1. Martin, K. E.
(2015). Synchrophasor measurements under the IEEE standard c37.118.1-2011 with amendment c37.118.1a. IEEE Transactions on Power Delivery, 30(3), 1514–1522. https://doi.org/10.1109/tpwrd.2015.2403591
↵
1. Mo, Y.,
2. Deng, X.,
3. Zheng, L., &
4. Liu, Z.
(2013). Kalman-consensus filter for time synchronization in wireless sensor networks. IET International Conference on Information and Communications Technologies (IETICT 2013). https://doi.org/10.1049/cp.2013.0080
↵
1. Myrda, P. T., &
2. Koellner, K.
(2010). NASPInet - the internet for synchrophasors. 2010 43rd Hawaii International Conference on System Sciences. https://doi.org/10.1109/hicss.2010.283
↵
1. Najafi, M.,
2. Nadealian, Z.,
3. Rahmanian, S., &
4. Ghafarinia, V.
(2019). An adaptive distributed approach for the real-time vision-based navigation system. Measurement, 145, 14–21. https://doi.org/10.1016/j.measurement.2019.05.015
↵
1. Nuñez, J. R.,
2. Anderton, C. R., &
3. Renslow, R. S.
(2018). Optimizing colormaps with consideration for color vision deficiency to enable accurate interpretation of scientific data ( J. Malo, Ed.). PLOS ONE, 13(7), e0199239. https://doi.org/10.1371/journal.pone.0199239
CrossRef
↵
1. Ochieng, W. Y.,
2. Sauer, K.,
3. Walsh, D.,
4. Brodin, G.,
5. Griffin, S., &
6. Denney, M.
(2003). GPS integrity and potential impact on aviation safety. Journal of Navigation, 56(1), 51–65. https://doi.org/10.1017/s0373463302002096
↵
1. Oldham, K. B.,
2. Myland, J. C., &
3. Spanier, J.
(2008). The error function erf(x) and its complement erfc(x). In An atlas of functions (pp. 405–415). Springer US. https://doi.org/10.1007/978-0-387-48807-3_41
↵
1. Paul, K.
(2011). Soft GNSS.
↵
1. Risbud, P.,
2. Gatsis, N., &
3. Taha, A.
(2018). Vulnerability analysis of smart grids to GPS spoofing. 2018 IEEE Power & Energy Society General Meeting (PESGM). https://doi.org/10.1109/pesgm.2018.8586646
↵
1. Schmidt, D.,
2. Radke, K.,
3. Camtepe, S.,
4. Foo, E., &
5. Ren, M.
(2016). A survey and analysis of the GNSS spoofing threat and countermeasures. ACM Computing Surveys, 48(4), 1–31. https://doi.org/10.1145/2897166
↵
1. Shepard, D. P.,
2. Humphreys, T. E., &
3. Fansler, A. A.
(2012). Evaluation of the vulnerability of phasor measurement units to GPS spoofing attacks. International Journal of Critical Infrastructure Protection, 5(3–4), 146–153. https://doi.org/10.1016/j.ijcip.2012.09.003
↵
1. Shi, D.,
2. Chen, T., &
3. Shi, L.
(2015). On set-valued kalman filtering and its application to event-based state estimation. IEEE Transactions on Automatic Control, 60(5), 1275–1290. https://doi.org/10.1109/tac.2014.2370472
↵
1. Shiryaev, V., &
2. Podivilova, E.
(2015). Set-valued estimation of linear dynamical system state when disturbance is decomposed as a system of functions. Procedia Engineering, 129, 252–258. https://doi.org/10.1016/j.proeng.2015.12.045
↵
1. Simon, D.
(2007). Optimal state estimation: Kalman, h [infinity], and nonlinear approaches. Choice Reviews Online, 44(06), 44–3334–44–3334. https://doi.org/10.5860/choice.44-3334
↵
1. Talebi, S. P., &
2. Werner, S.
(2018). Distributed kalman filtering: Consensus, diffusion, and mixed. 2018 IEEE Conference on Control Technology and Applications (CCTA). https://doi.org/10.1109/ccta.2018.8511492
↵
1. Tholomier, D.,
2. Kang, H., &
3. Cvorovic, B.
(2009). Phasor measurement units: Functionality and applications. 2009 Power Systems Conference. https://doi.org/10.1109/psamp.2009.5262468
↵
1. Wang, L.,
2. Fernandez, J.,
3. Burgett, J.,
4. Conners, R. W., &
5. Liu, Y.
(2008). An evaluation of network time protocol for clock synchronization in wide area measurements. 2008 IEEE Power and Energy Society General Meeting - Conversion and Delivery of Electrical Energy in the 21st Century. https://doi.org/10.1109/pes.2008.4596234
↵
1. Xu, H.,
2. Angrisano, A.,
3. Gaglione, S., &
4. Hsu, L.-T.
(2020). Machine learning based LOS/NLOS classifier and robust estimator for GNSS shadow matching. Satellite Navigation, 1(1). https://doi.org/10.1186/s43020-020-00016-w
↵
1. Zhan, L.,
2. Liu, Y.,
3. Yao, W.,
4. Zhao, J., &
5. Liu, Y.
(2017). Utilization of chip-scale atomic clock for synchrophasor measurements. 2017 IEEE Power & Energy Society General Meeting. https://doi.org/10.1109/pesgm.2017.8274514

In this issue

Download PDF

Article Alerts

Email Article

Citation Tools

Bookmark this article

Cited By...

No citing articles found.

Google Scholar

More in this TOC Section

Show more Original Article

Keywords

[1] ↵
Alanwar, A.,
Said, H., &
Althoff, M.
(2019). Distributed secure state estimation using diffusion kalman filters and reachability analysis. 2019 IEEE 58th Conference on Decision and Control (CDC). https://doi.org/10.1109/cdc40024.2019.9029929

[2] Alanwar, A.,

[3] Said, H., &

[4] Althoff, M.

[5] ↵
Al-Hammouri, A. T.,
Nordstrom, L.,
Chenine, M.,
Vanfretti, L.,
Honeth, N., &
Leelaruji, R.
(2012). Virtualization of synchronized phasor measurement units within real-time simulators for smart grid applications. 2012 IEEE Power and Energy Society General Meeting. https://doi.org/10.1109/pesgm.2012.6344949

[6] Al-Hammouri, A. T.,

[7] Nordstrom, L.,

[8] Chenine, M.,

[9] Vanfretti, L.,

[10] Honeth, N., &

[11] Leelaruji, R.

[12] ↵
Allnutt, J.,
Anand, D.,
Arnold, D.,
Goldstein, A.,
Li-Baboud, Y.-S.,
Martin, A.,
Nguyen, C.,
Noseworthy, R.,
Subramaniam, R., &
Weiss, M.
(2017). Timing challenges in the smart grid (tech. rep.). National Institute of Standards; Technology. https://doi.org/10.6028/nist.sp.1500-08

[13] Allnutt, J.,

[14] Anand, D.,

[15] Arnold, D.,

[16] Goldstein, A.,

[17] Li-Baboud, Y.-S.,

[18] Martin, A.,

[19] Nguyen, C.,

[20] Noseworthy, R.,

[21] Subramaniam, R., &

[22] Weiss, M.

[23] ↵
Althoff, M.
(2016). CORA 2016 manual. TU Munich, 85748.

[24] Althoff, M.

[25] ↵
Althoff, M., &
Dolan, J. M.
(2014). Online verification of automated road vehicles using reachability analysis. IEEE Transactions on Robotics, 30(4), 903–918. https://doi.org/10.1109/tro.2014.2312453

[26] Althoff, M., &

[27] Dolan, J. M.

[28] ↵
Althoff, M.,
Stursberg, O., &
Buss, M.
(2009). Safety assessment for stochastic linear systems using enclosing hulls of probability density functions. 2009 European Control Conference (ECC). https://doi.org/10.23919/ecc.2009.7074473

[29] Althoff, M.,

[30] Stursberg, O., &

[31] Buss, M.

[32] ↵
Bhamidipati, S., &
Gao, G. X.
(2018). Multiple GPS fault detection and isolation using a graph-SLAM framework. Proceedings of the 31st International Technical Meeting of The Satellite Division of the Institute of Navigation (ION GNSS+ 2018). https://doi.org/10.33012/2018.16030

[33] Bhamidipati, S., &

[34] Gao, G. X.

[35] ↵
Bhamidipati, S., &
Gao, G. X.
(2019). GPS multireceiver joint direct time estimation and spoofer localization. IEEE Transactions on Aerospace and Electronic Systems, 55(4), 1907–1919. https://doi.org/10.1109/taes.2018.2879532

[36] Bhamidipati, S., &

[37] Gao, G. X.

[38] ↵
Bhamidipati, S., &
Gao, G. X.
(2020). GPS spoofing mitigation and timing risk analysis in networked PMUs via stochastic reachability. Proceedings of the 33rd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2020). https://doi.org/10.33012/2020.17757

[39] Bhamidipati, S., &

[40] Gao, G. X.

[41] ↵
Bhamidipati, S.,
Iiyama, K.,
Mina, T., &
Gao, G.
(2022). Time-transfer from terrestrial gps for distributed lunar surface communication networks. IEEE Aerospace Conference.

[42] Bhamidipati, S.,

[43] Iiyama, K.,

[44] Mina, T., &

[45] Gao, G.

[46] ↵
Bhamidipati, S.,
Mina, T. Y., &
Gao, G. X.
(2018). GPS time authentication against spoofing via a network of receivers for power systems. 2018 IEEE/ION Position, Location and Navigation Symposium (PLANS). https://doi.org/10.1109/plans.2018.8373542

[47] Bhamidipati, S.,

[48] Mina, T. Y., &

[49] Gao, G. X.

[50] ↵
Bujorianu, L. M.
(2012). Applications of Stochastic Reachability. Stochastic Reachability Analysis of Hybrid Systems, 203–207. https://doi.org/10.1007/978-1-4471-2795-6_11

[51] Bujorianu, L. M.

[52] ↵
Coleman, M. J., &
Beard, R. L.
(2020). Autonomous clock ensemble algorithm for GNSS applications. NAVIGATION, 67(2), 333–346. https://doi.org/10.1002/navi.366

[53] Coleman, M. J., &

[54] Beard, R. L.

[55] ↵
Crassidis, J. L., &
Junkins, J. L.
(2011). Optimal estimation of dynamic systems. Chapman; Hall/CRC. https://doi.org/10.1201/b11154

[56] Crassidis, J. L., &

[57] Junkins, J. L.

[58] ↵
Dagle, J.
(2018). Importance of synchrophasor technology in managing the grid. In Power electronics and power systems (pp. 1–11). Springer International Publishing. https://doi.org/10.1007/978-3-319-89378-5_1

[59] Dagle, J.

[60] ↵
Ebinuma, T.
(2018). GPS-SDR-SIM.

[61] Ebinuma, T.

[62] ↵
Georgescu, A.,
Gheorghe, A. V.,
Piso, M.-I., &
Katina, P. F.
(2019). Critical space infrastructures. In Critical space infrastructures (pp. 21–36). Springer International Publishing. https://doi.org/10.1007/978-3-030-12604-9_2

[63] Georgescu, A.,

[64] Gheorghe, A. V.,

[65] Piso, M.-I., &

[66] Katina, P. F.

[67] ↵
Girard, A.
(2005). Reachability of uncertain linear systems using zonotopes. In Hybrid systems: Computation and control (pp. 291–305). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-540-31954-2_19

[68] Girard, A.

[69] ↵
Günther, C.
(2014). A survey of spoofing and counter-measures. Navigation, 61(3), 159–177. https://doi.org/10.1002/navi.65

[70] Günther, C.

[71] ↵
Heng, L.,
Work, D. B., &
Gao, G. X.
(2015). GPS signal authentication from cooperative peers. IEEE Transactions on Intelligent Transportation Systems, 16(4), 1794–1805. https://doi.org/10.1109/tits.2014.2372000

[72] Heng, L.,

[73] Work, D. B., &

[74] Gao, G. X.

[75] ↵
Hetet, S.
(2000). Signal-to-noise ratio effects on the quality of gps observations. Fredericton, New Brunswick: Department of Geodesy and Geomatics, University of New Brunswick.

[76] Hetet, S.

[77] ↵
Hofmann-Wellenhof, B.,
Lichtenegger, H., &
Collins, J.
(1992). Global positioning system. Springer Vienna. https://doi.org/10.1007/978-3-7091-5126-6

[78] Hofmann-Wellenhof, B.,

[79] Lichtenegger, H., &

[80] Collins, J.

[81] ↵
Jiang, Z., &
Groves, P. D.
(2012). NLOS GPS signal detection using a dual-polarisation antenna. GPS Solutions, 18(1), 15–26. https://doi.org/10.1007/s10291-012-0305-5

[82] Jiang, Z., &

[83] Groves, P. D.

[84] ↵
Khalajmehrabadi, A.,
Gatsis, N.,
Akopian, D., &
Taha, A. F.
(2018). Real-time rejection and mitigation of time synchronization attacks on the global positioning system. IEEE Transactions on Industrial Electronics, 65(8), 6425–6435. https://doi.org/10.1109/tie.2017.2787581

[85] Khalajmehrabadi, A.,

[86] Gatsis, N.,

[87] Akopian, D., &

[88] Taha, A. F.

[89] ↵
Khalife, J., &
Kassas, Z. M.
(2017). Modeling and analysis of sector clock bias mismatch for navigation with cellular signals. 2017 American Control Conference (ACC). https://doi.org/10.23919/acc.2017.7963500

[90] Khalife, J., &

[91] Kassas, Z. M.

[92] ↵
Kuusniemi, H.,
Lachapelle, G., &
Takala, J. H.
(2004). Position and velocity reliability testing in degraded GPS signal environments. GPS Solutions, 8(4), 226–237. https://doi.org/10.1007/s10291-004-0113-7

[93] Kuusniemi, H.,

[94] Lachapelle, G., &

[95] Takala, J. H.

[96] ↵
Li, W.,
Gong, D.,
Liu, M.,
Chen, J., &
Duan, D.
(2013). Adaptive robust kalman filter for relative navigation using global position system. IET Radar, Sonar & Navigation, 7(5), 471–479. https://doi.org/10.1049/iet-rsn.2012.0170

[97] Li, W.,

[98] Gong, D.,

[99] Liu, M.,

[100] Chen, J., &

[101] Duan, D.

[102] ↵
Martin, K. E.
(2015). Synchrophasor measurements under the IEEE standard c37.118.1-2011 with amendment c37.118.1a. IEEE Transactions on Power Delivery, 30(3), 1514–1522. https://doi.org/10.1109/tpwrd.2015.2403591

[103] Martin, K. E.

[104] ↵
Mo, Y.,
Deng, X.,
Zheng, L., &
Liu, Z.
(2013). Kalman-consensus filter for time synchronization in wireless sensor networks. IET International Conference on Information and Communications Technologies (IETICT 2013). https://doi.org/10.1049/cp.2013.0080

[105] Mo, Y.,

[106] Deng, X.,

[107] Zheng, L., &

[108] Liu, Z.

[109] ↵
Myrda, P. T., &
Koellner, K.
(2010). NASPInet - the internet for synchrophasors. 2010 43rd Hawaii International Conference on System Sciences. https://doi.org/10.1109/hicss.2010.283

[110] Myrda, P. T., &

[111] Koellner, K.

[112] ↵
Najafi, M.,
Nadealian, Z.,
Rahmanian, S., &
Ghafarinia, V.
(2019). An adaptive distributed approach for the real-time vision-based navigation system. Measurement, 145, 14–21. https://doi.org/10.1016/j.measurement.2019.05.015

[113] Najafi, M.,

[114] Nadealian, Z.,

[115] Rahmanian, S., &

[116] Ghafarinia, V.

[117] ↵
Nuñez, J. R.,
Anderton, C. R., &
Renslow, R. S.
(2018). Optimizing colormaps with consideration for color vision deficiency to enable accurate interpretation of scientific data ( J. Malo, Ed.). PLOS ONE, 13(7), e0199239. https://doi.org/10.1371/journal.pone.0199239
CrossRef

[118] Nuñez, J. R.,

[119] Anderton, C. R., &

[120] Renslow, R. S.

[121] ↵
Ochieng, W. Y.,
Sauer, K.,
Walsh, D.,
Brodin, G.,
Griffin, S., &
Denney, M.
(2003). GPS integrity and potential impact on aviation safety. Journal of Navigation, 56(1), 51–65. https://doi.org/10.1017/s0373463302002096

[122] Ochieng, W. Y.,

[123] Sauer, K.,

[124] Walsh, D.,

[125] Brodin, G.,

[126] Griffin, S., &

[127] Denney, M.

[128] ↵
Oldham, K. B.,
Myland, J. C., &
Spanier, J.
(2008). The error function erf(x) and its complement erfc(x). In An atlas of functions (pp. 405–415). Springer US. https://doi.org/10.1007/978-0-387-48807-3_41

[129] Oldham, K. B.,

[130] Myland, J. C., &

[131] Spanier, J.

[132] ↵
Paul, K.
(2011). Soft GNSS.

[133] Paul, K.

[134] ↵
Risbud, P.,
Gatsis, N., &
Taha, A.
(2018). Vulnerability analysis of smart grids to GPS spoofing. 2018 IEEE Power & Energy Society General Meeting (PESGM). https://doi.org/10.1109/pesgm.2018.8586646

[135] Risbud, P.,

[136] Gatsis, N., &

[137] Taha, A.

[138] ↵
Schmidt, D.,
Radke, K.,
Camtepe, S.,
Foo, E., &
Ren, M.
(2016). A survey and analysis of the GNSS spoofing threat and countermeasures. ACM Computing Surveys, 48(4), 1–31. https://doi.org/10.1145/2897166

[139] Schmidt, D.,

[140] Radke, K.,

[141] Camtepe, S.,

[142] Foo, E., &

[143] Ren, M.

[144] ↵
Shepard, D. P.,
Humphreys, T. E., &
Fansler, A. A.
(2012). Evaluation of the vulnerability of phasor measurement units to GPS spoofing attacks. International Journal of Critical Infrastructure Protection, 5(3–4), 146–153. https://doi.org/10.1016/j.ijcip.2012.09.003

[145] Shepard, D. P.,

[146] Humphreys, T. E., &

[147] Fansler, A. A.

[148] ↵
Shi, D.,
Chen, T., &
Shi, L.
(2015). On set-valued kalman filtering and its application to event-based state estimation. IEEE Transactions on Automatic Control, 60(5), 1275–1290. https://doi.org/10.1109/tac.2014.2370472

[149] Shi, D.,

[150] Chen, T., &

[151] Shi, L.

[152] ↵
Shiryaev, V., &
Podivilova, E.
(2015). Set-valued estimation of linear dynamical system state when disturbance is decomposed as a system of functions. Procedia Engineering, 129, 252–258. https://doi.org/10.1016/j.proeng.2015.12.045

[153] Shiryaev, V., &

[154] Podivilova, E.

[155] ↵
Simon, D.
(2007). Optimal state estimation: Kalman, h [infinity], and nonlinear approaches. Choice Reviews Online, 44(06), 44–3334–44–3334. https://doi.org/10.5860/choice.44-3334

[156] Simon, D.

[157] ↵
Talebi, S. P., &
Werner, S.
(2018). Distributed kalman filtering: Consensus, diffusion, and mixed. 2018 IEEE Conference on Control Technology and Applications (CCTA). https://doi.org/10.1109/ccta.2018.8511492

[158] Talebi, S. P., &

[159] Werner, S.

[160] ↵
Tholomier, D.,
Kang, H., &
Cvorovic, B.
(2009). Phasor measurement units: Functionality and applications. 2009 Power Systems Conference. https://doi.org/10.1109/psamp.2009.5262468

[161] Tholomier, D.,

[162] Kang, H., &

[163] Cvorovic, B.

[164] ↵
Wang, L.,
Fernandez, J.,
Burgett, J.,
Conners, R. W., &
Liu, Y.
(2008). An evaluation of network time protocol for clock synchronization in wide area measurements. 2008 IEEE Power and Energy Society General Meeting - Conversion and Delivery of Electrical Energy in the 21st Century. https://doi.org/10.1109/pes.2008.4596234

[165] Wang, L.,

[166] Fernandez, J.,

[167] Burgett, J.,

[168] Conners, R. W., &

[169] Liu, Y.

[170] ↵
Xu, H.,
Angrisano, A.,
Gaglione, S., &
Hsu, L.-T.
(2020). Machine learning based LOS/NLOS classifier and robust estimator for GNSS shadow matching. Satellite Navigation, 1(1). https://doi.org/10.1186/s43020-020-00016-w

[171] Xu, H.,

[172] Angrisano, A.,

[173] Gaglione, S., &

[174] Hsu, L.-T.

[175] ↵
Zhan, L.,
Liu, Y.,
Yao, W.,
Zhao, J., &
Liu, Y.
(2017). Utilization of chip-scale atomic clock for synchrophasor measurements. 2017 IEEE Power & Energy Society General Meeting. https://doi.org/10.1109/pesgm.2017.8274514

[176] Zhan, L.,

[177] Liu, Y.,

[178] Yao, W.,

[179] Zhao, J., &

[180] Liu, Y.

Main menu

User menu

Search

GPS Spoofing Mitigation and Timing Risk Analysis in Networked Phasor Measurement Units via Stochastic Reachability

Abstract

1 INTRODUCTION

1.1 Design Aspects of the Proposed SR-DKF Algorithm

1.1.1 A spatially-distributed network of static GPS receivers

1.1.2 A distributed processing framework

1.1.3 A set representation via p-Zonotopes

1.2 Contributions

2 PRELIMINARIES OF P-ZONOTOPES

3 SYSTEM MODEL

4 THE PROPOSED SR-DKF: A SET-VALUED SECURE STATE ESTIMATION ALGORITHM

4.1 During Initialization: Pre-Defined Process and Measurement Error Bounds in Non-Spoofed Conditions

4.2 Set-Valued DKF

4.2.1 Predicted p-Zonotope of the Local State Vector via Time Update

4.2.2 Corrected p-Zonotope of the State Estimation Error via Measurement Update

4.3 Spoofing attack mitigation and data transfer

4.4 Timing Risk Analysis

4.5 Implementation Details

5 RESULTS AND ANALYSIS

5.1 During a Coordinated Signal-Level Spoofing Attack

5.2 Robustness Analysis of the Proposed SR-DKF Algorithm

6 CONCLUSIONS

HOW TO CITE THIS ARTICLE

ACKNOWLEDGEMENTS

APPENDIX A PRELIMINARY FINDINGS OF THE POINT-VALUED DKF

REFERENCES

In this issue

Citation Manager Formats

Related Articles

Cited By...

More in this TOC Section

Similar Articles

Keywords

Main menu

User menu

Search

GPS Spoofing Mitigation and Timing Risk Analysis in Networked Phasor Measurement Units via Stochastic Reachability

Abstract

1 INTRODUCTION

1.1 Design Aspects of the Proposed SR-DKF Algorithm

1.1.1 A spatially-distributed network of static GPS receivers

1.1.2 A distributed processing framework

1.1.3 A set representation via p-Zonotopes

1.2 Contributions

2 PRELIMINARIES OF P-ZONOTOPES

3 SYSTEM MODEL

4 THE PROPOSED SR-DKF: A SET-VALUED SECURE STATE ESTIMATION ALGORITHM

4.1 During Initialization: Pre-Defined Process and Measurement Error Bounds in Non-Spoofed Conditions

4.2 Set-Valued DKF

4.2.1 Predicted p-Zonotope of the Local State Vector via Time Update

4.2.2 Corrected p-Zonotope of the State Estimation Error via Measurement Update

4.3 Spoofing attack mitigation and data transfer

4.4 Timing Risk Analysis

4.5 Implementation Details

5 RESULTS AND ANALYSIS

5.1 During a Coordinated Signal-Level Spoofing Attack

5.2 Robustness Analysis of the Proposed SR-DKF Algorithm

6 CONCLUSIONS

HOW TO CITE THIS ARTICLE

ACKNOWLEDGEMENTS

APPENDIX A PRELIMINARY FINDINGS OF THE POINT-VALUED DKF

REFERENCES

In this issue

Citation Manager Formats

Jump to section

Related Articles

Cited By...

More in this TOC Section

Similar Articles

Keywords