GPS Spoofing Mitigation and Timing Risk Analysis in Networked PMUs via Stochastic Reachability

The main contributions of this paper are listed as follows:

1. 1.

   Utilizing an efficient set representation known as probabilistic zonotope (p-Zonotope) [24] to parameterize the stochastic reachable set, we develop a set-valued DKF framework to compute secure GPS time for a network of receivers.

2. 2.

   Considering offline-estimated measurement error bounds in non-spoofed (or authentic) conditions, we design a two-tier approach that performs

   1. (a)

      Measurement-level spoofing attack mitigation by adaptively evaluating 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, such as meaconing, signal-level spoofing, coordinated attack, etc.

   2. (b)

      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 violating the IEEE C37.118.1a-2014 standards [9].

3. 3.

   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 improvement in the state estimation accuracy via the proposed SR-DKF algorithm as compared to conventional point-valued DKF and single-receiver-based adaptive KF. 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 magnitude of spoofing attacks.

The rest of the paper is organized as follows: Section 2 outlines the system model, the preliminaries of stochastic reachability and the conventional point-valued DKF; Section 3 describes the proposed SR-DKF algorithm; Section 4 experimentally validates the improved state estimation accuracy of the GPS timing computed via SR-DKF as compared to the point-valued DKF, and also demonstrates the trend of timing risk associated with the estimated GPS time; and Section 5 concludes the paper.

A p-Zonotope $\mathcal{L}$, as seen in Eq. (1a) and Fig. 3(a), is characterized by three parameters [24], namely, a) the mean of its center, denoted by $\bm{c}\in\mathbb{R}^{n}$, b) the uncertainty in the center, represented by the generator matrix $G\in\mathbb{R}^{n\times e}$, and c) the over-bounding covariance of the Gaussian tails, denoted by $\Sigma\in\mathbb{R}^{n\times n}$.

Similarly, a p-Zonotope with zero covariance $\Sigma=\bm{0}$ simplifies to a zonotope [26] that encompasses all the possible values of the center, such that the generator matrix of $Z$ is $G=[\bm{g}^{1},\cdots,\bm{g}^{e}]$. This case is mathematically represented by $Z$ in Eq. (1c).

The p-Zonotope is a combination of both sets, i.e., $\mathcal{L}=\mathcal{Z}\boxplus Z$, where $\boxplus$ is a set-addition operator, such that
$\mathcal{L}=sup\big{\{}f_{\mathcal{Z}}\Big{|}\mathbb{E}[\mathcal{Z}]\in Z\big{\}}$, where $f_{\mathcal{Z}}$ is the Probability Density Function (PDF) of the distributions represented by $\mathcal{Z}$ and $\mathbb{E}[\cdot]$ denotes the expectation operator. From Eqs. (1a)-(1c), note that $\mathcal{L}$ depends on both $G$ and $\underline{G}$. Additionally, note that unlike a PDF, the area enclosed by a p-Zonotope does not equal to one. More details regarding the characteristics of zonotopes and p-Zonotopes can be found in the prior literature [24, 23].

For example, the 2D p-Zonotope shown in Fig. 3(b) is formulated by considering the bounds on the mean and covariance of the first state to be $[-5,5]$ and $[0,2]$, respectively, and of the second state to be $[-10,10]$ and $[0,3]$, respectively. Mathematically, we represent this 2D p-Zonotope $\mathcal{L}=(\bm{c},G,\Sigma)$ by $\bm{c}=\begin{bmatrix}0\\ 0\end{bmatrix}$, $G=\begin{bmatrix}2&0\\ 0&3\end{bmatrix}$ and $\Sigma=3\times\begin{bmatrix}2&0\\ 0&3\end{bmatrix}=\begin{bmatrix}6&0\\ 0&9\end{bmatrix}$. We consider a multiplication factor $3$ in the over-bounding covariance $\Sigma$ based on the empirical rule according to which $99.7\%$ of the values lie within three standard deviations of the mean [27].

To perform set operations on the stochastic reachable sets, we require three important properties [24]:

• •

  Minkowski sum of addition: The addition of two p-Zonotopes $\mathcal{L}_{1}=(\bm{c}_{1},G_{1},\Sigma_{1})$ and $\mathcal{L}_{2}=(\bm{c}_{2},G_{2},\Sigma_{2})$ is given by

  	$\displaystyle\mathcal{L}_{1}\bigoplus\mathcal{L}_{2}$	$\displaystyle=(\bm{c}_{1}+\bm{c}_{2},[G_{1},G_{2}],\Sigma_{1}+\Sigma_{2}),$		(2a)
  	$\displaystyle\bigoplus_{i=1}^{p}\mathcal{L}_{i}$	$\displaystyle=(\sum_{i=1}^{p}\bm{c}_{i},[G_{1},\cdots,G_{p}],\sum_{i=1}^{p}\Sigma_{i}),$		(2b)
  where $\bigoplus$ is known as the Minkowski sum operator and $[\cdot,\cdot]$ denotes a horizontal concatenation operator.

• •

  Linear map: Given a matrix $A\in\mathbb{R}^{n\times n}$ and a p-Zonotope $\mathcal{L}=(\bm{c},G,\Sigma)$,

  $\displaystyle A\mathcal{L}$

  $\displaystyle=(A\bm{c},AG,A\Sigma A^{\top})$

  (2c)

• •

  Translation: Given a vector $\mu\in\mathbb{R}^{n}$ and a p-Zonotope $\mathcal{L}=(\bm{c},G,\Sigma)$,

  $\displaystyle\mu+\mathcal{L}$

  $\displaystyle=(\mu+\bm{c},G,\Sigma)$

  (2d)

Based on the key features of the SR-DKF algorithm, we mathematically represent the inter-connected network of power substations, each mounted with one GPS receiver, as a undirected graph $\mathcal{G}=\{M,E\}$, where $M$ indicates the set of nodes and $E$ denotes the set of connections between the nodes. Each node represents a GPS receiver in the network and each edge indicates the presence of a secure bidirectional communication link between the power substations. The total number of receivers in the network are $\big{|}M\big{|}$, where $\big{|}\cdot\big{|}$ denotes the cardinality of a set. Our SR-DKF algorithm does not require the network of power substations to be fully connected, and hence, $\big{|}E\big{|}\leq\dfrac{\big{|}M\big{|}\big{(}\big{|}M\big{|}-1\big{)}}{2}$.

The neighborhood set of any $i$th receiver $\forall i\in M$ is represented by $M_{i}$ and indicates the subset of nodes among the $M$ nodes that can communicate with the $i$th receiver, including itself. Therefore, the associated cardinality of the neighborhood set is denoted by $\big{|}M_{i}\big{|}$ where $\big{|}M_{i}\big{|}\geq 1\leavevmode\nobreak\ \forall i\in M$. Every receiver in the network has a prior knowledge regarding the pre-computed static position of all the receivers in its neighborhood set.

We define the global state vector at the $k$th time iteration, which is represented by $\bm{x}_{k}\in\mathbb{R}^{2}$, to be comprised of GPS time and its drift rate, which are given by $T_{k}$ and $\dot{T}_{k}$, respectively. As explained earlier in Section 1, estimated GPS time $T_{k}$ is used in PMUs for time-tagging the phasor measurements. At the $k$th iteration, the global state vector is the same across all the receivers in the network. Additionally, we define the local state vector at the $i$th receiver by $\bm{y}_{i,k}\in\mathbb{R}^{2}$. The local state vector comprises two states, namely, the receiver clock bias and clock drift, which are denoted by $c\delta t_{k}$ and $c\delta\dot{t}_{k}$, respectively. At the $i$th receiver, the global and local state vectors are related via secure time-varying local clock parameters, i.e., $T^{i}_{rx,k}$ and $\Delta T^{i}_{rx,k}$, which denote the measured local clock time and the sampling time, respectively. Therefore, $\bm{x}_{k}=\bm{y}_{i,k}+\bm{\theta}_{i,k}$ where $\bm{\theta}_{i,k}=[T^{i}_{rx,k},\Delta T^{i}_{rx,k}]^{\top}$. At the $k$th iteration, each $i$th receiver maintains an estimate of both the local and global state vectors. When a receiver is spoofed, the local state vector $\bm{y}_{i,k}$ is manipulated, thereby misleading the global GPS time $T_{k}$.

Based on [22], the general framework of a point-valued DKF implementation performed independently at the $i$th receiver, is shown in Eqs. (4)-(6). During initialization, i.e., $k<0$, we define the predicted value of the point-valued mean and covariance that are given by $\hat{\bm{x}}_{i,-1}$ and $\hat{P}_{i,-1}$, respectively. The process and measurement noise are modeled via a zero-mean Gaussian distribution that are given by $\nu_{i,k}\sim\mathcal{N}(\bm{0},Q_{i,k})$ and $\omega_{i,k}\sim\mathcal{N}(\bm{0},R_{i,k})$, 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: one is the measurement update and the other is the time update.

During the time update, as seen in Eq. (5), the corrected variables $\bar{\bm{x}}_{i,k-1}$ and $\bar{P}_{i,k-1}$ from the previous time, i.e., at the $(k-1)$th iteration, are propagated forward-in-time to predict the estimates for the next $k$th iteration, which are given by $\hat{\bm{x}}_{i,k}$ and $\hat{P}_{i,k}$, respectively.

At the $k$th iteration, each $i$th receiver receives the GPS measurements $\bm{z}_{j,k}$ and measurement noise covariance matrix $R_{j,k}$ from the other receivers in the neighborhood set, i.e., $j\in 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., $\bar{\bm{x}}_{i,k}$ and $\bar{P}_{i,k}$, respectively. The optimal Kalman gain denoted by $K_{j,k}\leavevmode\nobreak\ \forall j\in M_{i}$ is seen in Eq. (6a). For an adaptive implementation of point-valued DKF, additional parameters, namely forgetting factor $\psi$ and pre-measurement residuals $\bm{\epsilon}_{i,k}$, are considered in Eq. (6b) to adaptively estimate the measurement noise covariance matrix $R_{i,k}$ at each time iteration.

We formulate the set-valued DKF by applying stochastic reachability to the point-valued DKF described earlier in
Eqs. (4)-(5). Similar to point-valued DKF, the set-valued DKF also performs time and measurement updates to compute the predicted and corrected stochastic reachable sets of the global state vector $\bm{x}_{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. At the $k$th iteration, the estimation error in the point-valued predicted and corrected state estimate is given by $\Delta\hat{\bm{x}}_{i,k}$ and $\Delta\bar{\bm{x}}_{i,k}$, respectively, where $\Delta\hat{\bm{x}}_{i,k}=\hat{\bm{x}}_{i,k}-\bm{x}_{k}$ and $\Delta\bar{\bm{x}}_{i,k}=\bar{\bm{x}}_{i,k}-\bm{x}_{k}$.

At the $k$th iteration, we utilize the translation set property in Eq. (2d) to represent the predicted and corrected p-Zonotopes of the global state vector $\mathcal{L}^{i}_{\hat{\bm{x}},k}$ and $\mathcal{L}^{i}_{\bar{\bm{x}},k}$, respectively, as $\mathcal{L}^{i}_{\hat{\bm{x}},k}=\bm{x}_{k}+\mathcal{L}^{i}_{\Delta\hat{\bm{x}},k}$ and $\mathcal{L}^{i}_{\Delta\bar{\bm{x}},k}=\bm{x}_{k}+\mathcal{L}^{i}_{\Delta\bar{\bm{x}},k}$. Therefore, the predicted and corrected p-Zonotopes of the global state vector are given by Eq. (9). Note that translating the center of p-Zonotope as seen in Eq. (9) does not change the generator and covariance matrix. Therefore, the generator $G$ and covariance matrix $\Sigma$ of the predicted (or corrected) p-Zonotope of the state estimate $\bm{x}$ and state estimation error $\Delta\bm{x}$ are the same. From Eq. (9) and prior literature [16, 24], it is established that the center mean of predicted and corrected p-Zonotopes is equivalent to their point-valued DKF estimates.

Before measurement update at the $i$th receiver and $k$th iteration, we independently utilize the predicted p-Zonotope of the state estimation error $\mathcal{L}^{i}_{\Delta\bm{x},k}$ from Eq. (11) and the pre-computed measurement noise bound $\mathcal{L}^{i}_{\omega,k}$ to estimate the p-Zonotope of expected measurement innovation. We apply the set properties of p-Zonotopes in Eqs. (2a)-(2c) on the point-valued measurement innovation $\bm{\epsilon}_{i,k}=\bm{z}_{i,k}-H_{i}\hat{\bm{x}}_{i,k}$, which is defined earlier in Eq. (6). This p-Zonotope denoted by $\mathcal{L}^{i}_{\bm{\epsilon},k}$ is given by

The p-Zonotope of expected measurement innovation $\mathcal{L}^{i}_{\bm{\epsilon},k}$ provides a probability value corresponding to each point-valued measurement innovation $\bm{\epsilon}_{i,k}$ and this probability is denoted by $\alpha_{i,k}$. Based on this, for any $i$th receiver, we independently compute the attack status $\tilde{\alpha}_{i,k}$ of the received GPS measurements by normalizing the obtained probability value $\alpha_{i,k}$ with the probability value corresponding to the center mean of the p-Zonotope $\mathcal{L}^{i}_{\bm{\epsilon},k}$ and subtracting it from one. The measurement attack status $\tilde{\alpha}_{i,k}\in\mathbb{R}$ lies between $[0,1]$. As seen earlier in Eq. (7), the pre-computed process and measurement noise bounds are defined for non-spoofed conditions. Therefore, in a non-spoofed condition, the point-valued measurement innovation $\bm{\epsilon}_{i,k}$ and its p-Zonotope $\mathcal{L}^{i}_{\bm{\epsilon},k}$ of expected measurement innovation are in close agreement, and hence a low attack status $\tilde{\alpha}_{i,k}\approx 0$ is obtained. However, during spoofing, the point-valued measurement innovation does not comply with this estimated p-Zonotope, and therefore, a high value of $\tilde{\alpha}_{i,k}\approx 1$ is observed.

Thereafter, each $i$th receiver in the network $i\in M$, broadcasts its received GPS measurement $\bm{z}_{i,k}$ and the estimated attack status $\tilde{\alpha}_{i,k}$ to all the receivers in its neighborhood set $j\in M_{i}-i$. The set $M_{i}-i$ simply implies that the $i$th receiver need not broadcast GPS measurements to itself. Similarly, each $i$th receiver also receives the GPS measurements $\bm{z}_{j,k}$ and their associated attack status $\tilde{\alpha}_{i,k}$ from one or more receivers in its neighborhood set $j\in M_{i}-i$. In an intuitive sense, the measurement attack status obtained from each $j$th receiver in the neighborhood set indicates its belief in the received GPS measurement $\bm{z}_{j,k}$. Thereafter, the measurement attack status is used to compute the adaptive gain matrix $\tilde{K}_{j,k}$ as seen in Eq. (15). Given that the gain matrix is a value between $0$ and $1$, the Eq. (15) empirically represents a joint minimization of both spoofing probability and covariance of the point-valued state estimate. The formulated adaptive gain matrix $\tilde{K}_{j,k}$ is later utilized to perform the set-valued measurement update of SR-DKF, which was described earlier in Eq. (13).

We evaluate the probability of corrected p-Zonotope $\mathcal{L}^{i}_{\Delta\bar{x},k}$ of state estimation error derived in Eq. (13), to intersect an unsafe set, i.e., a set of states that violate the IEEE C37.118.1a-2014 standards. As shown in both Figs. 7(a) and 7(b), the unsafe set represents two gray strips, i.e., $\Delta\bar{T}_{i,k}\in[26.5\mu s,\infty]$ and $\Delta\bar{T}_{i,k}\in[-\infty,-26.5\mu s]$. Intuitively, this matches the definition of risk described earlier in Section 1, which denotes the probability of the state estimation error to exceed a safety AL.

To perform the timing risk analysis, we over-approximate the corrected p-Zonotope $\mathcal{L}^{i}_{\Delta\bar{x},k}$ by a finite series of polytopes [24]. An example of over-approximation of p-Zonotope by four polytopes $D_{1}$ to $D_{4}$ is shown in Fig. 7(a). Thereafter, we compute the intersection area of each $l$th polytope $D_{l}$ with the unsafe set $\mathcal{B}$. Figure 7(b) shows the top-view of over-approximated series of polytopes and an example intersection region of $D_{3}$ polytope with unsafe set $\mathcal{B}$.

The timing risk, seen in Eq. (16), is defined as the sum across polytopes $D_{l}$ of their intersection area, which is denoted by $\mathcal{A}(\cdot)$, with the unsafe set $\mathcal{B}$ times the max probability of $D_{l}$th polytope, which is denoted by $p_{max,l}$. The max probability of $D_{l}$th polytope is at its top face and equals the probability of corrected p-Zonotope at that level. For computational tractability, we represent the p-Zonotope by a finite number of polytopes and a $\gamma$-confidence set. The $\gamma$-confidence set [24] thresholds the tail of a p-Zonotope by a pre-defined parameter $\gamma$, such that $P[-\gamma<\mathcal{N}(0,1)<\gamma]=\textrm{erf}(\dfrac{\gamma}{\sqrt{2}})$, where erf denotes the error function [30]. From [24], the probability that the state estimation error lies outside the $\gamma$-confidence set of corrected p-Zonotope is $1-\textrm{erf}\Big{(}\dfrac{\gamma}{\sqrt{(}2)}\Big{)}^{2n}$.

Figure 8 summarizes the implementation steps of the proposed SR-DKF described in Sections 3.1-3.3. We assume that the $i$th receiver $\leavevmode\nobreak\ \forall i\in M$ is authentic during initialization. Therefore, at iteration $k=0$, we define the predicted p-Zonotope $\mathcal{L}^{i}_{\Delta\bar{\bm{x}},-1}$ of the state estimation error to be zero-mean center, with non-zero generator and over-bounding covariance matrices.

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. 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\leavevmode\nobreak\$s, the first spoofing attack manipulates the GPS time at $Rx:5$ to deviate at a rate of $100\leavevmode\nobreak\$ns/s, and the second spoofing attack, which occurs between $k=800-1300\leavevmode\nobreak\$s, deviates the GPS time of $Rx:1$ at a rate of $400\leavevmode\nobreak\$ns/s.

In Fig. 10(b), we show the performance of the point-valued state estimation via an adaptive DKF framework, where even the receivers that are not directly spoofed, e.g., $Rx:3$, exhibit high state estimation errors. Also, in Fig. 10(c), we show the state estimation accuracy achieved via a point-valued single-receiver-based adaptive KF. In both these point-valued approaches, the measurement noise covariance $R_{i,k}$ is adaptively estimated via Eq. (6) by setting the forgetting factor $\psi=0.3$. Based on the maximum error values listed in Table 1, we observe that both the approaches violate the IEEE C37.118.1a-2014 standards. In Fig. 10(d), we validate that our proposed SR-DKF algorithm provides a secure estimate of GPS timing with maximum error of $8.4\leavevmode\nobreak\ \mu$s across all receivers, and thereby complies with the IEEE standards at all times. The quantitative maximum error values of the receiver network for the above-described state estimation approaches are listed in Table 1.

In Fig. 11, we plot the attack status $\tilde{\alpha}_{i,k}$ of the GPS measurements $\bm{z}_{i,k}$ at each receiver in the network. In accordance with Section 3.2, we demonstrate that the proposed SR-DKF algorithm successfully detects and mitigates spoofing in the GPS measurements of $Rx:1$ and $5$ while accurately categorizing the other receivers, i.e., $Rx:2,3,4,6,7$ to be non-spoofed with attack status $\tilde{\alpha}_{i,k}\approx 0$ at all times. We observe that the spoofing attack with higher drift rate is mitigated quicker, i.e., $400$ ns/s in $Rx:1$, as compared to the one with lower drift rate, i.e., $100$ ns/s in $Rx:5$.

In Fig. 12, we showcase the timing risk associated with the estimated GPS timing, i.e., the probability of GPS time to violate the IEEE C37.118.1a-2014 standards. This is computed by analyzing the intersection of the corrected p-Zonotope $\mathcal{L}^{i}_{\Delta\bar{\bm{x}},k}$ with the unsafe set $\mathcal{B}$, as discussed in Section 3.3. An increase or jump in measure of timing risk is based on the adaptive gain matrix $\tilde{K}_{j,k}\leavevmode\nobreak\ \forall j\in M_{i}$ that in turn depends on the estimated attack status $\tilde{\alpha}_{j,k}$ of the receiver. Therefore, an increase in the timing risk for $Rx:5$ till $t=1040\leavevmode\nobreak\$s is consistent with the decrease in measurement reliability, which is evident from Fig 11. By comparing Fig. 10(a) and Fig. 12, we observe a direct correlation between the value of timing risk to the number of reliable measurements available for the calculation of corrected p-Zonotope. For instance, since none of the four neighbors of $Rx:4$ are spoofed, the estimated timing risk is consistently low, i.e., $\approx 10^{-7}$, for the entire experiment duration. In contrast, $Rx:3$ has significantly high timing risk $(\approx 0.3)$ because two out of its three neighbors are spoofed. Similarly, when only one out of three neighbors at $Rx:1$ is spoofed, the order of timing risk magnitude is $10^{-3}$ whereas when one out of two neighbors are spoofed at $Rx:5$, the timing risk is a higher order magnitude of $0.1$. Therefore, our SR-DKF estimates a robust measure of the timing risk across all receivers in the network.