Detection of False Data Injection Attack in Smart Grid via Adaptive Kalman Filtering
-
摘要: 研究了一种针对智能电网中假数据注入攻击的有效检测方法. 假数据注入攻击可以保持攻击前后残差基本不变, 绕过传统的不良数据检测技术. 首先基于电网模型, 分析了假数据注入攻击的攻击特性, 针对噪声统计特性未知且无迹Kalman滤波 (Unscented Kalman filter, UKF) 不稳定的现象, 提出了自适应平方根无迹Kalman滤波改进算法. 基于状态估计值, 结合中心极限定理提出检测算法, 并与欧几里得检测方法、巴氏系数检测方法进行比较. 最后, 仿真表明本文所提检测算法的优越性.
-
关键词:
- 智能电网 /
- 虚假数据注入攻击 /
- 攻击检测 /
- 自适应平方根无迹卡尔曼滤波
Abstract: In this paper, an effective detection method for false data injection attack in smart grid is studied. False data injection attack can keep the residual unchanged before and after the attack and bypass the traditional bad data detection technology. Firstly, based on the grid model, the attack characteristics of false data injection attack are analyzed. Aiming at the phenomenon that the noise statistical characteristics are unknown and the unscented Kalman filter (UKF) is unstable, an effective detection method for false data injection attack is proposed. An improved algorithm of adaptive square root unscented Kalman filter is proposed. Based on the state estimation and the central limit theorem, the algorithm is compared with Euclidean method and bayonet coefficient method. Finally, the simulation shows the superiority of the algorithm. -
智能电网是一种新型的电网, 它采用先进的通信网络技术和控制技术来支持更高效的能源安全传输和分配. 然而, 由于智能电网系统的复杂性和开放性, 智能电网中进行数据交换的信息网络成为易受到恶意攻击的对象[1-2]. 例如, 2016年, 黑客攻击乌克兰国家电力部门致使国内发生了一次大规模的停电事件[3], 造成严重经济损失. 因此, 智能电网的攻击检测研究具有重要意义.
隐蔽假数据攻击是目前恶意攻击的典型代表, 攻击者对传感器节点注入精心设计的错误数据, 接收错误数据的控制中心, 继而做出错误决策破坏系统的稳定性[4]. 在文献[5]中, 拒绝服务攻击旨在中断电力网络通信信道的可用性. 文献[6]构建了一种对智能电网中以完整性和可用性为目标的攻击分类方法. 文献[7]解决了在电力网络中攻击检测和拓扑隔离的问题. 文献[8]提出了一种算法来识别要操作的智能电表的最优数目, 从而找到最优攻击策略, 其目的是干扰电网系统的状态估计, 影响其稳定性. 文献[9]设计了具有隐蔽特性的虚假数据攻击, 它可以使攻击前后残差基本不变, 因此, 基于卡方检测器的检测技术是无效的. 近年来, 隐蔽性攻击检测成为了研究热点之一.
针对智能电网遭受虚假数据注入攻击的检测问题, 近年来有了很多成果. 文献[10]中, 以电机的电压模型为研究对象, 使用卡尔曼滤波(Kalman filter, KF)技术得到该系统的残差序列, 通过计算攻击前后残差序列之间的Bhattacharyya距离来判断系统中是否存在攻击. 该文献的不足之处有两点: 1) Bhattacharyya距离无法辨别虚假数据注入攻击前后两个残差序列的相似性, 因为虚假数据注入攻击前后残差保持不变; 2)该电压模型是线性的, 对于实际中存在的大多数的非线性系统, 该检测算法是失效的. 文献[11]使用卡尔曼滤波技术得到系统的状态, 从系统状态角度考虑设计检测算法, 提出欧几里德检测方法, 该方法可以检测隐蔽假数据注入攻击. 该文献的不足之处是, 系统方程对于噪声是线性的, 即该噪声是加性噪声, 当系统方程对于噪声是非线性的时候, 噪声经过非线性变换, 不再服从高斯分布, 也就无法设计阈值, 无法检测攻击. 文献[12]中, 攻击者可能会缓慢改变多个传感器的测量值, 因此上述统计异常检测不会检测到个别受损的测量值, 所以提出检测思想, 这些测量值组合起来会导致状态变量远离其真值, 然而, 文中作出假设, 控制中心收集到的测量值都是服从高斯分布的, 并基于此性质设计了双边假设检验检测攻击. 该文献的不足之处是, 电网系统中的参数仅会在特定情况下服从高斯分布, 而这种情况并不常见, 且建模时没有考虑噪声, 所以该文中的检测算法局限性较大.
因此, 根据虚假数据注入攻击的特性, 考虑系统非线性和噪声统计特性未知情况, 本文提出一种基于自适应平方根无迹卡尔曼滤波器 (Unscented Kalman filter, UKF) 的智能电网隐蔽假数据攻击检测方法. 该算法可以有效地、稳定地应用于非线性系统中对系统状态作出估计, 依据系统状态设计检测算法检测攻击, 并从检测指标的角度与现有算法进行对比. 最后进行仿真实验, 实验结果证明, 所提出的基于该算法的攻击检测方法可以准确地给出状态估计值, 从而检测出隐蔽攻击.
1. 电网模型、攻击特性及问题描述
1.1 电网模型
考虑一个由3条互连总线组成的电力网络, 与之对应的拓扑图
$G(V,\varepsilon )$ , 如图1所示.$V = \left\{ i \right\}_1^N$ 表示节点集,$i \in V$ 与总线$i$ 相对应,$\varepsilon \subseteq V \times V$ 是图$G$ 的边集. 对于总线$i \in V$ 的相角, 建立一个连续时间二阶微分动力学方程, 称为摆动方程[7]$${m_i}{\ddot \delta _i}(t) + {d_i}{\dot \delta _i}(t) = {u_i}(t) + {\eta _i}(t) - \sum\limits_{j \in {N_i}} {{P_{ij}}(t)} $$ (1) 式中,
${\delta _i}$ 表示总线的相角,${m_i}$ 为转动惯量,${d_i}$ 为阻尼系数,${u_i}$ 为控制输入信号,${\eta _i}$ 是过程噪声,${P_{ij}}$ 是从总线$i$ 到$j$ 的有功功率流$${P_{ij}}(t) = {V_i}{V_j}{b_{ij}}\sin ({\delta _i}(t) - {\delta _j}(t))$$ 式中,
$V$ 表示总线电压,${b_{ij}}$ 表示总线$i,j$ 之间的电纳.引入角速度
$\omega $ , 将上式改写为一阶微分方程组的形式, 即$$\begin{split} & {{\dot \delta }_i}(t) = {\omega _i}(t) \\ &{{\dot \omega }_i}(t) = - \frac{1}{{{m_i}}}({d_i}{\omega _i}(t) + {u_i} + {\varpi _i}(t) - \sum\limits_{j \in {N_i}} {{P_{ij}}(t)} ) \end{split} $$ (2) 每一条总线都安装有传感器, 将总线的数据传输回控制中心, 系统输出方程为
$${\boldsymbol{y}_{i}} = \left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right]{\boldsymbol{x}_{i}} + \boldsymbol{\rho }(t)$$ (3) 式中,
${\boldsymbol{y}_{i}}$ 表示总线${i}$ 的输出向量,$\boldsymbol{\rho }$ 是测量噪声.针对上述系统, 本文采用欧拉离散化方法[13], 即对于一个非线性系统
$$\dot{\boldsymbol x} = f(\boldsymbol{x}) + B\tau + \varpi $$ 离散化后的形式为
$${\boldsymbol{x}_{k}} = {\boldsymbol{x}_{{k} - 1}} + \tau f({\boldsymbol{x}_{{k} - 1}}) + \tau B\boldsymbol{u} + \tau \varpi + {\rm{O}}({\tau ^2})$$ 其中
$${\rm{O}}({\tau ^2}) = \dfrac{{{\tau ^2}}}{2}\dot f(\boldsymbol{x}) + \dfrac{{{\tau ^3}}}{3}{f^{(2)}}(\boldsymbol{x}) + \cdots + \dfrac{{{\tau ^{{{{k}}}}}}}{{k!}}{f^{({{{k}}} - 1)}}(\boldsymbol{x})$$ 采样间隔值
$\tau $ 取值充分小以保证高阶项${\rm{O}}({\tau ^2})$ 忽略不计, 本文中$\tau = 0.1\;{\rm{s}}$ .对式 (2) 使用上述欧拉离散化算法对该系统离散化, 得到
$$\left\{\begin{aligned} &{\boldsymbol{x}_{{k}}} = {f_{{o}}}({\boldsymbol{x}_{{k} - 1}}) + {B_o}\boldsymbol{u} + {\boldsymbol{\eta }_{k - 1}}\\\\ &{\boldsymbol{y}_{{k}}} = H{\boldsymbol{x}_{{k}}} + {\boldsymbol{\rho }_k}\\ \end{aligned} \right. $$ (4) 其中,
${f_o}( \cdot )$ 表示关于$x$ 的非线性函数, 将${\boldsymbol{x}_{{k} - 1}} + \tau {{f}}({\boldsymbol{x}_{{k} - 1}})$ 合并在一起,${B_o} = \tau B$ ,$\eta = \tau \varpi $ .假设
$\eta $ 和$\rho $ 是互不相关的非零均值高斯白噪声, 统计特性为$$\left\{ {\begin{array}{*{20}{l}} {E({{{\boldsymbol{\eta }}}_k}) = {{\boldsymbol{q}}},\;\;{{ }}{\rm{cov}} ({{{\boldsymbol{\eta }}}_k}{{\boldsymbol{\eta }}}_j^{\rm{T}}) = {Q_k}{\delta _{kj}}} \\ {E({{{\boldsymbol{\rho }}}_k}) = {{\boldsymbol{r}}},\;\;{{ }}{\rm{cov}} ({{{\boldsymbol{\rho }}}_k}{{\boldsymbol{\rho }}}_j^{\rm{T}}) = {R_k}{\delta _{kj}}} \\ {E({\eta _k}\rho _j^{\rm{T}}) = 0} \end{array}} \right.$$ (5) 式中, Q是非负定矩阵, R是正定矩阵,
${\delta _{kj}}$ 是$Krone cher{\text{-}} \delta$ 函数, 即$${\delta _{kj}} = \left\{ {\begin{aligned} &{0{{ }},\;k \ne j}\\ &{1{{ }},\;k = j} \end{aligned}} \right.$$ 将
$\eta $ 和$\rho $ 改写为${\eta _{k}} = \boldsymbol{q} + {\mu _{k}}$ ,${\rho _{{k}}} = \boldsymbol{r} + {\boldsymbol{v}_{{k}}}$ , 有$$\left\{ {\begin{array}{*{20}{l}} {{\boldsymbol{x}_{k}} = {f_o}({\boldsymbol{x}_{{k} - 1}}) + {B_o}\boldsymbol{u} + \boldsymbol{q} + {\boldsymbol{\rho }_{k - 1}}}\\ {{\boldsymbol{y}_{{k}}} = H{\boldsymbol{x}_{{k}}} + \boldsymbol{r} + {\boldsymbol{v}_{k - 1}}} \end{array}} \right.$$ (6) 式中,
${\mu _k}$ 和${{\boldsymbol{v}}_k}$ 是互不相关的零均值高斯白噪声.${\boldsymbol{x}_{{k}}} = {[{\delta _1} \quad {\delta _2}\quad {\delta _3}\quad {\omega _1}\quad {\omega _2}\quad {\omega _3}]^{\rm{T}}}$ ,${\boldsymbol{y}_{{k}}} \in {\bf{R}}^{6}$ .1.2 攻击描述和问题提出
假设一个恶意的第三方想要破坏第1.1节中描述的系统的完整性. 本文主要考虑虚假数据注入攻击.
假设攻击者具有系统知识, 知道系统矩阵, 控制矩阵, 测量矩阵, 可以控制一系列系统中的传感器数据.
考虑假数据注入攻击模型描述为
$${\boldsymbol{z}_{a}}(k) = \left\{ {\begin{aligned} &{H{\boldsymbol{x}_{{k}}} + \xi (k),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;k\not \subset T}\\ &{H{\boldsymbol{x}_{{a,k}}} + \xi (k) + {\boldsymbol{y}_{{a}}}(k),\;\;\;\;\;k \subset T} \end{aligned}} \right.$$ (7) 式中,
$\varGamma = {{\rm{diag}}}\{{{\gamma _1},\cdots,{\gamma _n}}\}$ 代表传感器选择矩阵, T是攻击的时间范围, 当${\gamma _i} = 1$ 时, 代表第$i$ 个传感器遭受攻击, 否则${\gamma _i} = 0$ . 且${\boldsymbol{y}_{{a}}}(k)$ 是攻击者精心构建的攻击向量序列, 攻击框图如图2所示.在隐蔽假数据注入攻击情况下, 在所考虑的攻击模型中, 攻击者构造攻击序列注入传感器, 在该序列下, 状态差值最终将发散到
$\infty $ , 而不会触发卡方检测器. 该攻击序列满足以下性质[9]:$$\begin{split} &\mathop {\lim }\limits_{k \to T} \left\| {{\boldsymbol{x}_{{a}}}(k) - \boldsymbol{x}(k)} \right\| > \alpha, \alpha > 0\\ &\mathop {\lim }\limits_{k \to T} \left\| {{\boldsymbol{r}_{{a}}}(k)} \right\| \le \varepsilon,\;\;\;\;\;\;\;\;\;\;\;\;\; 0 < \varepsilon < 1 \end{split}$$ (8) 其中,
${\boldsymbol{x}_{{a}}}({{k}})$ 和${\boldsymbol{r}_{{a}}}({{k}})$ 是遭受攻击时系统的状态变量和残差. 对攻击者克服检测机制的隐蔽性进行分析[4]. 由输出方程可得$$\begin{split} {\boldsymbol{r}_a} =& \left\| {(\boldsymbol{y} + \boldsymbol{a}) - H(\hat{\boldsymbol x} + \boldsymbol{c})} \right\| =\\ &\left\| {(\boldsymbol{y} - H\hat{\boldsymbol x}) + (\boldsymbol{a} - H\boldsymbol{c})} \right\|\leq \\ & \left\| {\boldsymbol{y} - H\hat{\boldsymbol x}} \right\| + \left\| {\boldsymbol{a} - H\boldsymbol{c}} \right\| = \boldsymbol{r} + {\tau _a} \end{split}$$ (9) 式中,
$\boldsymbol{a}$ 表示攻击信号,$\hat{\boldsymbol x}$ 表示状态估计值,$\boldsymbol{c}$ 表示状态变化量. 由式(9)知, 当$\boldsymbol{a} = H\boldsymbol{c}$ 时,${\boldsymbol{r}_{{a}}} = \boldsymbol{r}$ , 从输出残差角度考虑, 实现了隐蔽性.由上述分析可知, 从输出残差角度无法实现攻击检测. 所以, 需要借助滤波器从系统状态入手, 利用状态差值检测攻击信号. 本文考虑通过基于自适应平方根无迹卡尔曼滤波器的方法, 从系统状态入手来检测此种攻击.
2. 自适应平方根无迹卡尔曼滤波
本文提出将噪声估计环节[14]加入到无迹卡尔曼滤波算法中, 实现对于噪声统计特性的在线估计, 同时改变标准UKF中状态误差协方差矩阵的迭代方式来保证滤波器的稳定性.
2.1 平方根无迹变换
本小节给出平方根无迹变换的实现方式.
步骤 1. 利用
${{k}} - 1$ 时刻的估计状态及状态误差协方差矩阵来计算sigma采样点$\xi {}_i{{ }}$ , 并给出其权值${W_{m,i}}$ 和${W_{c,i}}$ .$$\left\{ {\begin{array}{*{20}{l}} {{\xi _{\boldsymbol{i},k - 1}} = {{\hat{\boldsymbol x}}_{k - 1}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1}\\ {{\xi _{\boldsymbol{i},k - 1}} = {{\hat{\boldsymbol x}}_{k - 1}} + \gamma {{({S_{k - 1}})}_i},\;\;\;i = 2, \cdots ,n + 1}\\ {{\xi _{\boldsymbol{i},k - 1}} = {{\hat{\boldsymbol x}}_{k - 1}} - \gamma {{({S_{k - 1}})}_i},\;\;\;i = n + 2, \cdots ,2n + 1} \end{array}} \right.$$ 式中,
$n$ 代表状态的维数,$\gamma = \sqrt {n + \lambda } $ . 这里不再借助${\hat P_{k - 1}}$ 构造采样点, 而用${S_{k - 1}}$ .${S_{k - 1}}$ 是对状态误差协方差矩阵$P$ 进行QR分解得到的$P$ 的平方根, QR分解返回的是$P$ 的Cholesky分解的平方根的转置, 即${(\sqrt {(n + \lambda )P} )^{\rm{T}}}$ , 故${({\boldsymbol{S}_{k - 1}})_i}$ 等于$(\sqrt {(n + \lambda )P} )_i^{\rm{T}}$ , 这两种sigma点的选取方法在数学上是等价的.$${W_{i,m}} = \left\{ {\begin{array}{*{20}{l}} {\frac{\lambda }{{n + \lambda }}{{ }},\;\;\;\;\;i = 1} \\\\ {\frac{1}{{2(n + \lambda )}},\;{{ }}i = 2, \cdots ,2n + 1} \end{array}} \right.\;\;\;\;\;\;\;\;$$ $$\;\;\;\;\;\;\;{W_{i,c}} = \left\{ {\begin{array}{*{20}{l}} {\frac{\lambda }{{n + \lambda }} + {{(1}} + \beta - {\alpha ^2}{{) }},i = 1} \\\\ {\frac{1}{{2(n + \lambda )}}{{ }},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 2, \cdots ,2n + 1} \end{array}} \right.$$ 式中,
$\lambda = {\alpha ^2}(n + \kappa ) - n$ 用来控制采样点之间的距离, 一般${10^{ - 4}} \leq \alpha \leq 1$ ,$\beta = 2$ ,$n + \kappa = 3$ [8].步骤 2. 计算sigma点通过非线性函数的传播结果
$$\begin{split} &{\tau _i} = f({\xi _i}),\;\;\;i = 1, \cdots ,2n + 1\\ &\boldsymbol{z} = \sum\limits_{i = 1}^{2n + 1} {{W_{m,i}}{\boldsymbol{\tau }_i}} \\ &{P_x} = \sum\limits_{i = 1}^{2n + 1} {{W_{c,i}}({\boldsymbol{\tau }_i} - \boldsymbol{z}){{({\boldsymbol{\tau }_i} - \boldsymbol{z})}^{\rm{T}}}} \end{split}$$ (10) 式(10)是无迹变换中的计算误差协方差矩阵的公式, 平方根无迹变换将式(10)换为
$$\begin{split}{S_z} =\;& \boldsymbol{qr}(\sqrt {{\boldsymbol{W}_{c,1}}} ({\tau _1} - \boldsymbol{z})\times \cdots\times \\ &\sqrt {{\boldsymbol{W}_{c,2{n} + 1}}} ({\tau _{2n + 1}} - \boldsymbol{z}))\end{split}$$ (11) 2.2 噪声估计方法
假设噪声是非零均值的高斯噪声, 本文使用的是Sage-Husa噪声估计器, 通过极大后验估计原理, 获得次优噪声估计值, 噪声估计部分可参考文献[14].
对于系统(6), 系统方程是非线性的, 输出方程是线性的, 所以噪声估计器为
$${\hat{\boldsymbol q}_{k + 1}} = \frac{1}{k}\left[(k - 1){\hat{\boldsymbol q}_k} + {\hat{\boldsymbol x}_{k + 1}} - \sum\limits_{i = 1}^{2n + 1} {W_i^m} {f_o}({\xi _{i,k}})\right]\quad$$ $$\begin{split} {{\hat Q}_{k + 1}} =& \dfrac{1}{k}\big[({k} - 1){{\hat Q}_k} + {K_{k + 1}}{\boldsymbol{\varepsilon }_{k + 1}}{\boldsymbol{\varepsilon }_{k + 1}}K_{k + 1}^{\rm{T}} + {P_{k + 1}}-\\ & \sum\limits_{i = 1}^{2n + 1} {W_i^c} ({\boldsymbol{\tau }_{i,k + 1\left| k \right.}} - {{\hat{\boldsymbol x}}_{k + 1\left| k \right.}})\times\\ &{({\boldsymbol{\tau }_{i,k + 1\left| k \right.}} - {{\hat{\boldsymbol x}}_{k + 1\left| k \right.}})^{\rm{T}}}\big] \end{split}$$ $$\begin{split} &{{\hat{\boldsymbol r}}_{k + 1}} = \dfrac{1}{k}\left[(k - 1){{\hat{\boldsymbol r}}_k} + {\boldsymbol{y}_{k + 1}} - H{{\hat{\boldsymbol x}}_{k + 1\left| k \right.}}\right]\\ &{{\hat{\boldsymbol R}}_{k + 1}} = \dfrac{1}{k}\left[(k - 1){{\hat{\boldsymbol R}}_k} + {\boldsymbol{\varepsilon }_{k + 1}}\boldsymbol{\varepsilon }_{k + 1}^{\rm{T}} - H{P_{k + 1\left| k \right.}}H_k^{\rm{T}}\right] \end{split}$$ 式中,
${\boldsymbol{\varepsilon }_k} = {\boldsymbol{y}_k} - {\hat{\boldsymbol y}_{k\left| {k - 1} \right.}}$ 代表残差,$K$ 代表滤波器增益矩阵.2.3 自适应平方根无迹卡尔曼滤波设计
本节用到了MATLAB中的3个函数, 分别为: qr函数, cholupdate函数, diag函数.
1)
$qr\;( \cdot )$ 表示QR分解,$A \in {\mathbb{{\boldsymbol{{\rm{R}}}}}^{m \times l}}$ ,$[Q,R] = $ $qr(A)$ , 该函数生成$m \times l$ 上三角形矩阵R和$m \times m$ 酉矩阵Q, 从而,$A = Q\times R$ .2)
$R{{1}} = cholupdate(R,x)$ , 返回$A + xx'$ 的上三角Cholesky因子,$x$ 是具有合适长度的一个列向量, 其中,$R = chol(A)$ 是A的原始Cholesky分解因子.3)
${\rm{diag}}\{ \cdot \}$ 函数,$D = {\rm{diag}}\{v\}$ , 返回以向量v的元素为主对角线的对角矩阵.$x ={\rm{ diag}}\{A\}$ , 返回A的主对角线元素的列向量. 对于式(22)中最外侧的${\rm{diag}} \{\cdot \}$ , 返回一个对角矩阵, 主对角线元素是每一次迭代得到的噪声矩阵主对角线元素的平方根的实时估计值.自适应平方根无迹卡尔曼滤波(Adaptive square-root UKF, ASRUKF)算法的具体步骤如下.
步骤 1.初始化
$(k\;=\;1):$ $$\left\{ {\begin{array}{*{20}{l}} {{{\hat{\boldsymbol x}}_1} = {\rm{E}}({\boldsymbol{x}_1})}\\ {{S_1} = chol\left\{ {\rm{E}}\left(({\boldsymbol{x}_1} - {{\hat{\boldsymbol x}}_1}){{({\boldsymbol{x}_1} - {{\hat{\boldsymbol x}}_1})}^{\rm{T}}}\right)\right\} }\\ {\sqrt {{{\hat Q}_1}} = {S_1}}\\ {{{\hat{\boldsymbol q}}_1} = 0} \end{array}} \right.$$ (12) 步骤 2. 对于
$k\;=\;2,3,\cdots,$ 进行迭代a)时间更新
计算sigma点, 构造矩阵
$${\boldsymbol{\xi }_{k - 1}} = \left[{\hat{\boldsymbol x}_{k - 1}}{\hat{\boldsymbol x}_{k - 1}} + \gamma {({S_{k - 1}})_i}{\hat{\boldsymbol x}_{k - 1}} - \gamma {({S_{k - 1}})_i}\right]$$ (13) $$\left\{ {\begin{array}{*{20}{l}} {{\tau _{i,k\left| {k - 1} \right.}} = f({\boldsymbol{\xi }_{i,k - 1}}) + \boldsymbol{q}}\\ {\boldsymbol{x} = \sum\limits_{i = 1}^{2n + 1} {{W_{m,i}}{\tau _{i,k\left| {k - 1} \right.}} + \hat{\boldsymbol q}} }\\ {{S_{k\left| {k - 1} \right.}} = qr\left\{ \left[ \sqrt {{W_{c,i}}} ({\tau _{2:2n + 1}} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}})\sqrt {{{\hat Q}_{k - 1}}} \right]\right\} }\\ {S = cholupdate}\left\{ \left[{S_{k\left| {k - 1} \right.}},\right.\right. \\ \;\;\;\;\;\;\left.\left.\sqrt {{W_{c,i}}} ({\tau _{2:2n + 1}} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}), \pm \right]\right\} \end{array}} \right.$$ (14) b)测量更新
$${\tilde{\boldsymbol y}_k} = {\boldsymbol{y}_k} - H{\hat{\boldsymbol x}_{k\left| {k - 1} \right.}}$$ (15) 计算混合误差协方差矩阵平方根
$${P_{xy}} = S_{k\left| {k - 1} \right.}^{\rm{T}}{S_{k\left| {k - 1} \right.}}{H_k}$$ (16) 计算新息协方差矩阵
$${S_{y,k}} = qr{\bigg\{ \left[H{S_{k\left| {k - 1} \right.}}){{ }}\sqrt R \right]^{\rm{T}}}\bigg\} $$ (17) 计算滤波增益
$${K_k} = \frac{\frac{P_{xy,k}} {S_{y,k}^{\rm{T}}}}{S_{y,k}}$$ (18) 更新状态值
$${\hat{\boldsymbol x}_k} = {\hat{\boldsymbol x}_{k\left| {k - 1} \right.}} + {K_k}{\tilde{\boldsymbol y}_k}$$ (19) 更新状态误差协方差矩阵平方根
$$\left\{ {\begin{array}{*{20}{l}} {U = {K_k}S_{y,k}^{\rm{T}}} \\\\ {{S_k} = cholupdate\left\{ {S_{k\left| {k - 1} \right.}}{{ }}U{{ }} - 1\right\} } \end{array}} \right.$$ (20) 估计过程噪声均值
$${\hat{\boldsymbol q}_{k + 1}} = \dfrac{1}{k}\left[(k - 1){\hat{\boldsymbol q}_k} + {\hat{\boldsymbol x}_{k + 1}} - \sum\limits_{i = 1}^{2n + 1} {W_i^m} {f_o}({\boldsymbol{\xi }_{i,k}})\right]$$ (21) 估计过程噪声协方差矩阵平方根
$$\left\{ {\begin{array}{*{20}{l}} {\sqrt {Q1} = cholupdate\left\{ \sqrt {{{\hat Q}_{k - 1}}} ,\left| {{{\hat{\boldsymbol x}}_k} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}} \right|,\dfrac{1}{k}\right\} }\\ {\sqrt {Q2} = cholupdate\left\{ \sqrt {Q1} ,U, - \dfrac{1}{k}\right\} }\\ {\sqrt {{{\hat Q}_k}} = {\rm{diag}}\left\{ \sqrt {{\rm{diag}}\left\{\sqrt {Q2} {{\sqrt {Q2} }^{\rm{T}}}\right\}} \right\} }\\[-10pt] \end{array}} \right.$$ (22) 注1. Sage-Husa噪声估计器不能同时处理过程噪声和测量噪声都未知的情况, 否则会造成滤波发散, 故上述算法假设测量噪声的均值和方差均是已知的[15].
注2. 式(14)的最后一个子式中,
$ \pm $ 表示当${W_{c,1}}> 0$ 时, 取正; 反之, 取负, 以保证根号下不出现负数[16].2.4 自适应平方根无迹卡尔曼滤波稳定性
本小节给出自适应平方根无迹卡尔曼滤波算法的估计误差保持随机有界性的充分条件, 并给出证明.
对于非线性电网系统(6)
$$\begin{split} &{\boldsymbol{x}_k} = {f_o}({\boldsymbol{x}_{k - 1}}) + {B_o}\boldsymbol{u} + \boldsymbol{q} + {\boldsymbol{\mu }_{k - 1}}\\ &{\boldsymbol{y}_k} = H{\boldsymbol{x}_{k - 1}} + \boldsymbol{r} + {\boldsymbol{v}_k} \end{split}$$ (23) 定义
$${\tilde{\boldsymbol x}_k} = {\boldsymbol{x}_k} - {\hat{\boldsymbol x}_k}\;\;\;\;\;\;\;\;\;\;\;\;$$ (24) $${\tilde{\boldsymbol x}_{k\left| {k - 1} \right.}} = {\boldsymbol{x}_k} - {\hat{\boldsymbol x}_{k\left| {k - 1} \right.}}$$ (25) 在
${\hat{\boldsymbol x}_{k - 1}}$ 处对${{{\boldsymbol{x}}}_k}$ 泰勒级数展开$$\begin{split}{\boldsymbol{x}_k} = & {\;f_o}({{\hat{\boldsymbol x}}_{k - 1}}) + {B_o}\boldsymbol{u} + \boldsymbol{q} + \nabla {f_o}({{\hat{\boldsymbol x}}_{k - 1}}){{\tilde{\boldsymbol x}}_{k - 1}}+\\ &\dfrac{1}{2}{\nabla ^2}{f_o}({{\hat{\boldsymbol x}}_{k - 1}})\tilde{\boldsymbol x}_{k - 1}^2 + \cdots + \omega \end{split}$$ (26) 其中
$${\nabla ^i}{f_{{o}}}(\hat{\boldsymbol x}){\tilde{\boldsymbol x}^i} = {\left(\sum\limits_{j = 1}^n {{{\tilde{\boldsymbol x}}_j}} \frac{\partial }{{\partial {\boldsymbol{x}_j}}}\right)^i}{\left. {{f_{{o}}}(\boldsymbol{x})} \right|_{\boldsymbol{x} = {{\hat{\boldsymbol x}}_{k - 1}}}}$$ 在
${\hat x_{k - 1}}$ 处对${\hat x_{k\left| {k - 1} \right.}}$ 泰勒级数展开$$\begin{split} {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}} = &\left(1 - \frac{n}{{n + \lambda }}\right) \times (f({{\hat{\boldsymbol x}}_{k - 1}}) + B\boldsymbol{u}) + \boldsymbol{q}\; + \\ &\frac{1}{{2(n + \lambda )}} \times \sum\limits_{i = 2}^{n + 1} {(f\left[({{\hat{\boldsymbol x}}_{k - 1}}) + \gamma {{({S_{k - 1}})}_i}\right] + } \\ &B\boldsymbol{u} + \boldsymbol{q}) + \frac{1}{{2(n + \lambda )}} \times \sum\limits_{i = 2}^{n + 1} {\bigg(f\bigg[({{\hat{\boldsymbol x}}_{k - 1}})\;- } \\ &\gamma {({S_{k - 1}})_i}\bigg] + B\boldsymbol{u} + \boldsymbol{q}\bigg)=\\ \;& f({{\hat{\boldsymbol x}}_{k - 1}}) + B\boldsymbol{u} + \boldsymbol{q} + {\nabla ^2}f({{\hat{\boldsymbol x}}_{k - 1}}){P_{k - 1}}\\[-10pt] \end{split}$$ (27) 将式(26)和式(27)代入式(25), 得
$${\tilde{\boldsymbol x}_{k\left| {k - 1} \right.}} = {\beta _k}{F_k}{\tilde{\boldsymbol x}_{k - 1}} + {\boldsymbol{\mu }_k}$$ (28) 式中
$${F_k} = {\left. {\dfrac{{\partial {f_o}({\boldsymbol{x}_{k - 1}})}}{{\partial {\boldsymbol{x}_{k - 1}}}}} \right|_{\boldsymbol{x} = {{\hat{\boldsymbol x}}_{k - 1}}}}$$ (29) 其中,
${\beta _k} = {\rm{diag}}\left\{ {{\beta _{1k}}, \cdots ,{\beta _{nk}}} \right\}$ 为未知的矩阵, 其作用是弥补建模一阶线性化引起的模型误差.对于式(24)中的输出方程, 定义
$${\tilde{\boldsymbol y}_{k\left| {k - 1} \right.}} = {\boldsymbol{y}_k} - {\hat{\boldsymbol y}_{k\left| {k - 1} \right.}}$$ (30) 式(30)可写成
$$\begin{split} {{\tilde{\boldsymbol y}}_{k\left| {k - 1} \right.}} =\;& H{\boldsymbol{x}_k} + \boldsymbol{r} + {\boldsymbol{v}_k} - (H{{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}} + \boldsymbol{r}) = \\ &H{{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} + \boldsymbol{r} + {\boldsymbol{v}_k} \end{split}$$ (31) 由式(14)可知,
${S_{k\left| {k - 1} \right.}}$ 是${\hat P_{k\left| {k - 1} \right.}}$ 的平方根, 根据$qr( \cdot )$ 函数的模式$[Q,R] = qr(A)$ , 此时$$A = {\left[\sqrt {{W_{c,i}}} \left({\boldsymbol{\tau }_{2:2n + 1,k\left| {k - 1} \right.}} - {\hat{\boldsymbol x}_{k\left| {k - 1} \right.}}\right)\sqrt {{{\hat Q}_{k - 1}}} \right]^{\rm{T}}}$$ 使用
$qr( \cdot )$ 函数有$A = Q\times R$ , 那么$$\begin{split} {R^{\rm{T}}}R =\; & R{Q^{\rm{T}}}QR = {A^{\rm{T}}}A =\\ & \sum\limits_{i = 2}^{2n + 1} {\omega _i}{{\left({\boldsymbol{\tau }_{i,k\left| {k - 1} \right.}} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}\right)}^{\rm{T}}}\; \times \\ &\left({\boldsymbol{\tau }_{i,k\left| {k - 1} \right.}} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}\right) + {{\hat Q}_k} \end{split}$$ 再根据
$cholupdate$ 函数的模式, 此时$$\begin{split} S_{k\left| {k - 1} \right.}^{\rm{T}}S_{k\left| {k - 1} \right.}^{} = \;\;& {R^{\rm{T}}}R + \sqrt {{W_{c,1}}} \left({\boldsymbol{\tau }_{1,k\left| {k - 1} \right.}} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}\right) \times \\ &\sqrt {{W_{c,1}}} {\left({\boldsymbol{\tau }_{1,k\left| {k - 1} \right.}} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}\right)^{\rm{T}}}=\\ &P_{k\left| {k - 1} \right.}^{}\\[-10pt] \end{split}$$ (32) 与之对比, 先验状态误差协方差矩阵的实际值为
$$\begin{split} {P_{k\left| {k - 1} \right.}} = \;&{\rm{E}}\left({{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}}\tilde{\boldsymbol x}_{k\left| {k - 1} \right.}^{\rm{T}}\right)=\\ & {\rm{E}}\left(({\beta _k}{F_k}{{\hat{\boldsymbol x}}_{k - 1}} + {\boldsymbol{\mu }_k}){({\beta _k}{F_k}{{\hat{\boldsymbol x}}_{k - 1}} + {\boldsymbol{\mu }_k})^{\rm{T}}}\right)=\\ & {\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + \Delta {P_{k\left| {k - 1} \right.}} + {Q_{k - 1}} \end{split}$$ (33) 式中,
$\Delta {P_{k\left| {k - 1} \right.}}$ 是${\rm{E}}\left(({\beta _k}{F_k}{\tilde{\boldsymbol x}_{k - 1}}){({\beta _k}{F_k}{\tilde{\boldsymbol x}_{k - 1}})^{\rm{T}}}\right)$ 与${\beta _k}\times {F_k} {\hat P_{k - 1}} F_k^{\rm{T}}{\beta _k}$ 的差值. 引入$\delta {P_{k\left| {k - 1} \right.}}$ 来表示误差协方差矩阵的真实值${P_{k\left| {k - 1} \right.}}$ 与${\hat P_{k\left| {k - 1} \right.}}$ 的差异, 式(32)可写为$$\begin{split} {{\hat P}_{k\left| {k - 1} \right.}} =\;& {P_{k\left| {k - 1} \right.}} + \delta {P_{k\left| {k - 1} \right.}} = \\ &{\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + \Delta {P_{k\left| {k - 1} \right.}} + \\ &Q{}_k +\delta {P_{k\left| {k - 1} \right.}}=\\ &{\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {\Xi _{k - 1}} \\ \end{split} $$ (34) 此时,
$\delta {P_{k\left| {k - 1} \right.}}$ 中包含了${Q_k}$ 与${\hat Q_k}$ 的差值, 因为在ASRUKF算法中噪声信息是未知的. 其中,${\Xi _k} = \Delta {P_{k\left| {k - 1} \right.}} + {Q_k} + \delta {P_{k\left| {k - 1} \right.}}$ .预测测量协方差矩阵
${\hat P_{yy}}$ 、互协方差矩阵${\hat P_{xy}}$ 以及后验状态误差协方差矩阵${\hat P_k}$ 可以写成$${\hat P_{yy}} = H{\hat P_{k\left| {k - 1} \right.}}{H^{\rm{T}}} + {R_k}\;\;\;\;$$ (35) $${\hat P_{xy}} = {\hat P_{k\left| {k - 1} \right.}}{H^{\rm{T}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$ (36) $${\hat P_k} = {\hat P_{k\left| {k - 1} \right.}} - {\hat P_{xy}}\hat P_{yy}^{ - 1}\hat P_{xy}^{\rm{T}}$$ (37) 引理 1[17]. 对于一个随机变量
${\xi _k}$ 及其函数$V({\xi _k})$ , 实数${\upsilon _{\min }},{\upsilon _{{\rm{\max}} }},\mu > 0,\;0 < \lambda \leq 1$ 对于任意的k, 若满足以下两个不等式性质:$$\begin{split} &{\upsilon _{\min }}{\left\| {{\boldsymbol{\xi }_k}} \right\|^2} \leq V({\xi _k}) \leq {\upsilon _{\max }}{\left\| {{\boldsymbol{\xi }_k}} \right\|^2}\\ &{\rm{E}}(V({\boldsymbol{\xi }_k})\left| {{\boldsymbol{\xi }_{k - 1}}} \right.) - V\left({\boldsymbol{\xi }_{k - 1}}\right) \leq \mu - \lambda V\left({\boldsymbol{\xi }_{k - 1}}\right) \end{split}$$ (38) 则该随机变量均方有界, 即
$$\begin{split} {\rm{E}}\left({\left\| {{\xi _k}} \right\|^2}\right) \leq& \dfrac{{{\upsilon _{\max }}}}{{{\upsilon _{\min }}}}{\rm{E}}\left({\left\| {{\xi _0}} \right\|^2}\right){(1 - \lambda )^k} + \\ &\dfrac{\mu }{{{\upsilon _{\min }}}}\sum\limits_{i = 1}^{k - 1} {{{(1 - \lambda )}^i}} \\ \end{split} $$ (39) 引理 2[17]. 假设矩阵
$A \in {{\boldsymbol{{\rm{R}}}}^{n \times m}}$ ,$B \in {{\boldsymbol{{\rm{R}}}}^{m \times n}}$ ,$C \in {{\boldsymbol{{\rm{R}}}}^{n \times n}}$ , 如果$A > 0,\;C > 0$ , 则有$${A^{ - 1}} > B{({B^{\rm{T}}}AB + C)^{ - 1}}{B^{\rm{T}}}$$ (40) 引理 3[17]. 假设矩阵
$A \in {{\boldsymbol{{\rm{R}}}}^{n \times n}}$ ,$C \in {{\boldsymbol{{\rm{R}}}}^{n \times n}}$ , 如果$A > 0$ ,$C > 0$ , 则有$${A^{ - 1}} > {(A + C)^{ - 1}}$$ (41) 定理 1. 对于非线性系统(23), 如果有如下假设条件成立
1)对于任意的k, 存在实数
${f_{\min }},{h_{\min }},{\beta _{\min }},$ ${f_{\max }}, {h_{\max }},{\beta _{\max }} \ne {{0}}$ , 使得下式成立:$$\begin{split} f_{\min }^2I \leq\; &{F_k}F_k^{\rm{T}} \leq f_{\max }^2I \\ h_{\min }^2I \leq\;& {H_k}H_k^{\rm{T}} \leq h_{\max }^2I \\ \beta _{\min }^2I \leq \;&{\beta _k}\beta _k^{\rm{T}} \leq \beta _{\max }^2I \\ \end{split} $$ (42) 2)存在实数
${\Xi _{\max }},{\Xi _{\min }},{r_{\min }},{p_{\max }},{p_{\min }}$ , 使得下式成立:$$\begin{split} & {p_{\min }}I \leq {P_k} \leq {p_{{\rm{max}}}}I \\ & {\Xi _{\min }}I < {\Xi _k} \leq {\Xi _{\max }}I \\ &{r_{\min }}I \leq {R_k} \\ \end{split} $$ (43) 那么自适应平方根无迹卡尔曼滤波算法的状态估计误差将是均方有界的, 即自适应平方根无迹卡尔曼滤波器稳定.
证明. 选择函数
$${V_k}({\tilde x_k}) = \tilde x_k^{\rm{T}}\hat P_k^{ - 1}{\tilde x_k}$$ (44) 由式(43)可知
$$\dfrac{1}{{{p_{\max }}}}{\left\| {{{\tilde x}_k}} \right\|^2} \leq V({\tilde x_k}) \leq \dfrac{1}{{{p_{\min }}}}{\left\| {{{\tilde x}_k}} \right\|^2}$$ 满足了引理1中第1个性质, 为满足引理1, 还需要得到式(38)中第2个性质
${\rm{E}}[{V_k}({x_k})] - {V_{k - 1}}({\tilde x_{k - 1}})$ 的上界.后验状态误差协方差矩阵
${\hat P_k}$ 的另一种表达式为$${\hat P_k} = {\hat P_{k\left| {k - 1} \right.}} - {\hat P_{xy}}\hat P_{yy}^{ - 1}\hat P_{xy}^{\rm{T}} = (I - {K_k}{H_k}){\hat P_{k\left| {k - 1} \right.}}$$ (45) 其中
$${K_k} = {\hat P_{k\left| {k - 1} \right.}}H_k^{\rm{T}}{({H_k}{\hat P_{k\left| {k - 1} \right.}}H_k^{\rm{T}} + {R_k})^{ - 1}}\;\;\;\;\;\;\;\;\;\;\;\;$$ (46) $${\tilde x_k} = {x_k} - ({\hat x_{k\left| {k - 1} \right.}} + {\hat P_{xy}}\hat P_{yy}^{ - 1}{\tilde y_k}) = {\tilde x_{k\left| {k - 1} \right.}} - {K_k}{\tilde y_k}\;\;$$ (47) 将式(47)代入式(44), 得
$$\begin{split} {V_k}({{\tilde{\boldsymbol x}}_k}) =\;& {\left({{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} - {K_k}{{\tilde{\boldsymbol y}}_k}\right)^{\rm{T}}}\hat P_k^{ - 1}\left({{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} - {K_k}{{\tilde{\boldsymbol y}}_k}\right) =\\ & {{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}}\hat P_k^{ - 1}{{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} - {\left[{H_k}{{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} + {v_k}\right]^{\rm{T}}}\; \times \\ &K_k^{\rm{T}}\hat P_k^{ - 1}{{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} - \tilde x_{k\left| {k - 1} \right.}^{\rm{T}}\hat P_k^{ - 1}{K_k}\big[{H_k}{{\tilde x}_{k\left| {k - 1} \right.}} + \\ &{v_k}\big] + {\left[{H_k}{{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} + {v_k}\right]^{\rm{T}}}K_k^{\rm{T}} \times \\ &\hat P_k^{ - 1}{K_k}\left[{H_k}{{\tilde{\boldsymbol x}}_{k\left| {k - 1} \right.}} + {v_k}\right]\\[-10pt] \end{split}$$ (48) 式(46)的另一种等价形式为
$${K_k} = (I - {K_k}{H_k}){\hat P_{k\left| {k - 1} \right.}}H_k^{\rm{T}}R_k^{ - 1} = {\hat P_k}H_k^{\rm{T}}R_k^{ - 1}$$ (49) 式(45)用矩阵求逆公式得
$$\hat P_k^{ - 1} = \hat P_{k\left| {k - 1} \right.}^{ - 1} + H_k^{\rm{T}}R_k^{ - 1}{H_k}$$ (50) 将式(49), 式(50)和式(28)代入式(48)中, 式(48)的条件期望为
$$\begin{split} {\rm{E}}({V_k}({{\tilde x}_k})\left| {{{\tilde x}_{k - 1}}} \right.) =\;& {\rm{E}}{\bigg(} {({A_{k - 1}}{{\tilde x}_{k - 1}} + {\mu _k})^{\rm{T}}} \times\\ &\hat P_{k\left| {k - 1} \right.}^{ - 1}({A_{k - 1}}{{\tilde x}_{k - 1}} + {\mu _k})-\;\;\\ & {\left[{B_k}\left({A_{k - 1}}{{\tilde x}_{k - 1}} + {\mu _k}\right)\right]^{\rm{T}}}\\ &\left(\Sigma _k^{ - 1} - \Sigma _k^{ - 1}{B_k}{{\hat P}_k}B_k^{\rm{T}}\Sigma _k^{ - 1}\right) \times \\ &\left[{B_k}({A_{k - 1}}{{\tilde x}_{k - 1}} + {\mu _k})\right] + \\ &v_k^{\rm{T}}\Sigma _k^{ - 1}{B_k}{{\hat P}_k}B_k^{\rm{T}}\Sigma _k^{ - 1}{v_k}\left| {{{\tilde x}_{k - 1}}} \right.{\bigg)} \end{split}$$ (51) 将式(34)代入式(51)中, 并由引理3得
$$\begin{split} &{{\rm{E}}{\bigg(}{V_k}({{\tilde{\boldsymbol x}}_k})\left| {{{\tilde{\boldsymbol x}}_{k - 1}}} \right. {\bigg)} } \leq\\ &\;\;\;{ {\rm{E}}{\bigg(} {{{({\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})}^{\rm{T}}}} {{\left( {{\beta _k}{F_k}\hat P_{k{\rm{ - }}1}^{}F_k^{\rm{T}}{\beta _k}} \right)}^{ - 1}}}\; \times \\ &\;\;\;{({\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}}) + \omega _k^{\rm{T}}{{({\beta _k}{F_k}\hat P_{k{\rm{ - }}1}^{}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})}^{ - 1}} \times }\\ &\;\;\;{{\omega _k} - {{({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})}^{\rm{T}}} \times \left[ {{H_k}} \right.\left( {{\beta _k}} \right.{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}} \times }\\ &\;\;\;{{\beta _k} + \left. {{{\hat Q}_k}} \right)H_k^{\rm{T}} + {{\left. {{R_k}} \right]}^{ - 1}}({H_k}{\beta _k}{F_k} \times {{\tilde{\boldsymbol x}}_{k - 1}})\;- }\\ &\;\;\;{{{({H_k}{\omega _k})}^{\rm{T}}}\left[ {{H_k}} \right.\left({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k}\right) \times H_k^{\rm{T}} + }\\ &\;\;\;{{{\left. {{R_k}} \right]}^{ - 1}}({H_k}{\omega _k}) + v_k^{\rm{T}}R_k^{ - 1}{H_k}{{\hat P}_k}H_k^{\rm{T}}R_k^{ - 1}{v_k} {\left| {{{\tilde{\boldsymbol x}}_{k - 1}}} \right.{\bigg)} }} \end{split}$$ (52) 观察式(52)的第1项
$$\begin{split} &{\rm{E}}{\bigg(} {({\beta _k}{F_k}{{\tilde x}_{k - 1}})^{\rm{T}}}{\left( {{\beta _k}{F_k}\hat P_{k{\rm{ - }}1}^{}F_k^{\rm{T}}{\beta _k}} \right)^{ - 1}} \times \\ &\qquad({\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})\left| {{{\tilde{\boldsymbol x}}_{k - 1}}} \right.{\bigg)} = \\ & \qquad\tilde{\boldsymbol x}_{k - 1}^{\rm{T}}{{\hat P}_{k - 1}}{{\tilde{\boldsymbol x}}_{k - 1}} = {V_{k - 1}}({{\tilde{\boldsymbol x}}_{k - 1}}) \end{split}$$ (53) 两边同时减去式(53), 得
$$\begin{split} &{{\rm{E}}\bigg({V_k}({{\tilde{\boldsymbol x}}_k})\left| {{{\tilde{\boldsymbol x}}_{k - 1}}} \right.\bigg)- {V_{k - 1}}({{\tilde {\boldsymbol{x}}}_{k - 1}})}\leq\\ &\;\;\;{ {\rm{E}}\bigg( {\omega _k^{\rm{T}}} {{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})}^{ - 1}} - }\\ &\;\;\;{{{({H_k}{\omega _k})}^{\rm{T}}} \times \left[ {{H_k}} \right.\big( {{\beta _k}} {F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} \;+ }\\ &\;\;\;{ {{{\hat Q}_k}} \big) H_k^{\rm{T}} + {{\left. {{R_k}} \right]}^{ - 1}}({H_k}{\omega _k}) + v_k^{\rm{T}}R_k^{ - 1}{H_k} \;\times }\\ &\;\;\;{{{\hat P}_k}H_k^{\rm{T}}R_k^{ - 1}{v_k} {\left| {{{\tilde{\boldsymbol{ x}}}_{k - 1}}} \right.} \bigg) - {{({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})}^{\rm{T}}} \times }\\ &\;\;\;{\left[ {{H_k}} \right.({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})H_k^{\rm{T}}}+\\ &\;\;\;{ {{\left. {{R_k}} \right]}^{ - 1}} \times \left({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}}\right)} \end{split}$$ (54) 观察式(54)的最后1项, 根据引理3, 得
$$\begin{split} &{\hat P_{k - 1}^{ - 1} > {{({H_k}{\beta _k}{F_k})}^{\rm{T}}}\left[ {({H_k}{\beta _k}{F_k}){{\hat P}_{k - 1}}{{({H_k}{\beta _k}{F_k})}^{\rm{T}}}} \right.}+\\ &\;\;{{{\left. { {H_k}{{\hat Q}_k}H_k^{\rm{T}} + {R_k}} \right]}^{ - 1}}({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})}\\[-10pt] \end{split}$$ (55) 式(55)分别左右乘
${\tilde {\boldsymbol{x}}_{k - 1}}$ , 有$$\begin{split} &\tilde{\boldsymbol x}_{k - 1}^{\rm{T}}\hat P_{k - 1}^{ - 1}{{\tilde {\boldsymbol{x}}}_{k - 1}} > {({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})^{\rm{T}}}[({H_k}{\beta _k}{F_k}){{\hat P}_{k - 1}} \times\\ &\quad{({H_k}{\beta _k}{F_k})^{\rm{T}}} + {H_k}{\Xi _k}H_k^{\rm{T}} + {R_k}{]^{ - 1}}({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}}) \end{split}$$ (56) 选择
$$\begin{split} {\lambda _k} =& \; \Bigl\{ {({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})^{\rm T}}\Bigl[({H_k}{\beta _k}{F_k}){{\hat P}_{k - 1}}\; \times \\ &{({H_k}{\beta _k}{F_k})^{ {\rm{T}}}} + {H_k}{Q _k}H_k^{\rm T} + {R_k}{\Bigl]^{ - 1}}\;\times\\ & ({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})\Bigl\}/(\tilde{\boldsymbol x}_{k - 1}^{\rm{T}}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}}) \end{split}$$ (57) 由式(56)和式(57)可知,
${\lambda _k} < 1$ , 得$$\begin{split} {\lambda _k} \geq \;&{p_{\min }}{({h_{\min }}{\beta _{\min }}{f_{\min }})^2} \; \times \\ &\big[{p_{\max }}{({h_{\max }}{\beta _{\max }}{f_{\max }})^2} + \\ &{{\hat q}_{\max }}h_{\max }^2 + {r_{\max }}{\big]^{ - 1}} = {\lambda _{\min }} > 0 \end{split} $$ (58) 由式(57), 得
$$\begin{split} &- {({H_k}{\beta _k}{F_k})^{\rm{T}}}[({H_k}{\beta _k}{F_k}){{\hat P}_{k - 1}}{({H_k}{\beta _k}{F_k})^{\rm{T}}}+\\ &\quad\quad\quad{H_k}{{\hat Q}_k}H_k^{\rm{T}} + {R_k}{]^{ - 1}}({H_k}{\beta _k}{F_k}{{\tilde{\boldsymbol x}}_{k - 1}})\leq \\ &\quad\quad\quad-{\lambda _{\min }}{V_{k - 1}}({{\tilde{\boldsymbol x}}_{k - 1}}) \end{split}$$ (59) 考虑式(54)中的其他项, 得
$$\begin{split} &{{\mu _k} = {\rm{E}}\bigg({\omega _k^{\rm{T}}}\Big[{{{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})}^{ - 1}}} }-\\ & \quad \; \quad{ {H_k}^{\rm{T}}{\left( {{H_k}} \right.}({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})H_k^{\rm{T}} + }\\ &\quad \; \quad{{\left. {{R_k}} \right)^{ - 1}} \times H_k\Big]}{\omega _k} + \\ &\quad \; \quad{\upsilon _k^{\rm{T}}R_k^{ - 1}{H_k}{{\hat P}_k}H_k^{\rm{T}}R_k^{ - 1}{v_k} {\left| {{{\tilde x}_{k - 1}}} \right.} \bigg)} \end{split}$$ (60) 考虑式(60)两边都是标量, 取迹不会改变等式, 所以根据
$${\rm{tr}}(AB) = {\rm{tr}}(BA)$$ (61) 得到
$$\begin{split} &{{\mu _k} = {\rm{tr}}\bigg\{ {\Big[ {{{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})}^{ - 1}}} } - {H^{\rm{T}}_k} \times }\\ &\quad\quad\;{\left( {{H_k}} \right.({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})H_k^{\rm{T}} + }\\ &\quad\quad\;{{{ {{R_k}} )}^{ - 1}}{ {{H_k}} \Big]} {{Q_k}} \bigg\} + {\rm{tr}}\{[R_k^{ - 1}{H_k}{{\hat P}_k}H_k^{\rm{T}}]\left| {{{\tilde {\boldsymbol{x}}}_{k - 1}}} \right.\} } \end{split}$$ (62) 根据引理2, 得
$$\begin{split} &{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})^{ - 1}} - \\ &\;\;\;\quad {H^{\rm{T}}_k}[{H_k}({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + \\ &\;\;\;\quad {{\hat Q}_k})H_k^{\rm{T}} + {R_k}{)^{ - 1}}{H_k}] > 0 \end{split}$$ (63) 由式(63)可知
${\mu _k} > 0$ , 由式(42)和式(43), 得$$\begin{split} {\mu _k} \leq\; & {\rm{tr}}[{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^{\rm{T}}{\beta _k} + {{\hat Q}_k})^{ - 1}}{{\hat Q}_k}]+\\ & {\rm{tr}}(R_k^{ - 1}{H_k}{{\hat P}_k}H_k^{\rm{T}})\leq\\ & \dfrac{{{q_{\max }}}}{{{{\hat q}_{\min }}}} \times L + \dfrac{{{h_{{{\max }^2}}}{p_{\max }}}}{{{r_{\min }}}} \times M = {\mu _{\max }} \end{split}$$ (64) 因此根据引理1, 得
$${\rm{E}}({V_k}({x_k})) - {V_{k - 1}}({\tilde x_{k - 1}}) \leq {\mu _{\max }} - {\lambda _{\min }}{V_{k - 1}}({\tilde x_{k - 1}})$$ (65) 保证是
${\tilde x_k}$ 有界的, 即$$\begin{split} {\rm{E}}({\left\| {{{\tilde{\boldsymbol x}}_k}} \right\|^2}) \leq\; & \dfrac{{{p_{\max }}}}{{{p_{\min }}}}{\rm{E}}({\left\| {{{\tilde{\boldsymbol x}}_0}} \right\|^2}){(1 - {\lambda _{\min }})^k}\;+\\ & \dfrac{{{\mu _{\max }}}}{{{p_{\min }}}}\sum\limits_{i = 1}^{k - 1} {{{(1 - {\lambda _{\min }})}^i}} \end{split}$$ (66) 因此, 自适应平方根无迹卡尔曼滤波算法状态估计误差将是均方有界的, 即自适应平方根无迹卡尔曼滤波器是稳定的. □
注3. 因为误差序列
$\left\{ {{{\hat{\boldsymbol x}}_k},k \geq 1} \right\}$ 是随机过程,${\rm{E}}({\left\| {{{\tilde{\boldsymbol x}}_k}} \right\|^2})$ 是该随机过程的均值函数, 随机有界指该随机过程有界. 即对于任意时刻$k$ , 对于服从高斯分布的噪声$\omega $ , 定理1确保了误差${\tilde{\boldsymbol x}_k}$ 的2范数平方的期望值是有界的, 即${\rm{E}}({\left\| {{{\tilde{\boldsymbol x}}_k}} \right\|^2}) \leq \alpha $ 上限$\alpha $ 的值由定理1确定.3. 检测方法
本节分析了文献[6], [10-11]中所使用的检测方法的不足之处. 提出一种新的检测方法, 利用中心极限定理设计新的检测算法, 依据随机变量的统计特性提出双边假设检验法检测攻击, 从误检率方面对本文算法给出分析, 并与巴氏系数检测算法[10]、欧几里德检测算法[11]进行对比与分析.
文献[6]针对虚假数据注入攻击, 构造了CUSUM算法来检测攻击. 该算法步骤如下:
$${\chi _t} = ({y_t} - {H_k}{\hat x_{t\left| {t - 1} \right.}})\Sigma _t^{ - 1}({y_t} - {H_k}{\hat x_{t\left| {t - 1} \right.}})\;$$ $${p_t} = 1 - F({\chi _t})\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$ $${s_t} = \log \left(\dfrac{\alpha }{{{p_t}}}\right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$ $$\begin{split} & \Gamma = \inf \left\{ {t:{g_t} \geq h} \right\} \\ &{g_t} = \max \left\{ {0,{g_{t - 1}} + {s_t}} \right\} \\ \end{split}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; $$ 文献[6]中假设噪声服从高斯分布, 可知
${\chi _t}$ 服从${\chi ^2}$ 分布.$F( \cdot )$ 表示${\chi ^2}$ 分布的概率分布函数,${p_t}$ 表示${\chi ^2}$ 分布的上(右)尾概率,$\alpha $ 表示系统给定的上(右)尾概率用于检验${y_t}$ 是否是异常值. 根据${\chi ^2}$ 分布的上$\alpha $ 分位点性质, 当${y_t}$ 是异常值时,${s_t} > 0$ ; 当${y_t}$ 是正常值时,${s_t} \leq 0$ . 当${y_t}$ 是异常值时,${g_t}$ 的值会随时间累积直到超过检测阈值$h$ ($h$ 的确定参考文献[6]), 此时发出警报. 该算法通过引入$h$ , 避免了欧氏检测算法那种仅考虑误检率的阈值设计方法, 在误检率和误检周期之间做到了平衡, 既降低了噪声的影响, 又可以尽快检测出攻击.但是, 该算法实质上是依赖与噪声的高斯分布特性, 即, 实质上还是
${\chi ^2}$ 检测, 所以在遇到本文所研究的攻击时, 攻击前后的残差基本保持不变, 因此该算法是不可用的.文献[10]中针对具有隐蔽性的复数域攻击, 提出一种检测思想, 滤波器在不受攻击下工作和滤波器在错误数据注入攻击下工作会产生两个服从高斯分布的残差序列, 通过计算上述两个服从高斯分布的残差序列的巴氏相似性系数
$DB$ 来判断攻击是否发生.当物理系统的基本动力学用线性系统和高斯统计特性的噪声进行建模时, 在正常操作期间, 残差序列
${v^o}$ 遵循零均值高斯分布${v^o} = {\rm{N}}({r^o},{S^o})$ , 攻击发生后残差序列$v_k^a = {\rm{N}}(r_k^a,S_k^a)$ .$$\begin{split} DB({v^o},v_k^a) =\;& \dfrac{1}{4}[R\{ {({r^o},r_k^a)^H}{\Lambda _k}({r^o},r_k^a)\}\;+ \\ & R\{ {({r^o},r_k^a)^H}{{\tilde \Lambda }_k}({r^o},r_k^a)\} ]\;+\\ & \dfrac{1}{2}\ln \frac{{\left| {{\Gamma _k}} \right|}}{{\sqrt {\left| {{S^o}} \right|\left| {S_k^a} \right|} }} \\ \end{split}$$ $${\Gamma _k} = \frac{1}{2}({S^o} + S_k^a) = \left[ {\begin{array}{*{20}{c}} {{\Gamma _k}}&{{{\tilde \Gamma }_k}} \\ {\tilde \Gamma _k^*}&{\Gamma _k^*} \end{array}} \right]\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$ 该检测方法的优点同时也是缺点, 即巴氏相似性系数
$DB$ 能检测复数域攻击, 但也仅能检测复数域攻击. 因为系统中有无复数域攻击, 差异主要体现在$DB({v^o},v_k^a)$ 的后两项上. 而本文的虚假数据注入攻击属于实数域范畴, 在正常操作期间, 残差序列${v^o}$ 遵循零均值高斯分布${v^o} = {\rm{N}}({m_1},{P_1})$ , 攻击发生后残差序列$v_k^a = {\rm{N}}({m_2},{P_2})$ , 实数域的巴氏距离计算式为$$\begin{split} DB(v_k^o,v_k^a) = \;&\dfrac{1}{8}[{({m_1} - {m_2})^{\rm{T}}}{P^{ - 1}}({m_1} - {m_2}) ]\;+\\ & \dfrac{1}{2}\ln \dfrac{{\det P}}{{\sqrt {\det {P_1}\times \det {P_2}} }} \\ \end{split} $$ $$P = \frac{{{P_1} + {P_1}}}{2}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$ Bhattacharyya距离无法辨别虚假数据注入攻击前后两个残差序列的相似性, 因为虚假数据注入攻击前后残差保持不变. 除此之外, 该检测法的前提是线性系统, 因为线性运算虽然会改变系统中随机变量的均值和方差但并不会改变随机变量服从高斯分布的特性, 对于本文所研究的非线性系统, 这种方法是不可用的.
目前存在的检测方法大多是在加性噪声的基础上提出的, 在文献[11]中, 假设过程方程和量测方程相对于噪声是线性的, 即
$$\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{x}_{k+1}} = f({\boldsymbol{x}_k},{u_k},{t_k}) + {\omega _k}$$ $$\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{y}_{k+1}} = h({\boldsymbol{x}_k},{t_k}) + {v_k}$$ 文献[11]中, 欧氏检测法的阈值和误检率正是基于噪声项的附加性得到的, 因为噪声并不参与到非线性变换中, 所以噪声的特性就是假设的正态分布的特性, 使用正态分布的
$3\sigma $ 原则设计阈值. 只有在这种前提下, 欧氏距离检测算法才是正确的. 实际上, 过程和量测方程相对于噪声也可能是非线性的, 即$$\qquad\qquad{\boldsymbol{x}_k} = f({\boldsymbol{x}_k},{u_k},{t_k},{\omega _k})$$ $$\qquad\qquad y_{k}=h(x_{k},t_{k},v_{k})$$ 在这种情况下, 欧氏距离检测法的阈值设计思想是行不通的, 因为噪声不再作为附加项, 而是直接参与到非线性变换中. 此时, 噪声不再服从正态分布, 统计特性未知, 因此阈值无法设计.
综上, 本文提出一种检测思想: 在攻击未发生时, 利用本文提出的非线性滤波方法, 根据总线
$i$ 上传感器在时间范围$[{t_0},{t_k}]$ 收集到的测量数据获得系统状态估计值, 进而获得后验状态误差值${\tilde{\boldsymbol x}_k}$ , 此时并不知道状态误差${\tilde{\boldsymbol x}_k}$ 这个随机变量服从什么分布, 而且也不用对该随机变量的统计特性作出假设. 本文将时间区间$[{t_0},{t_k}]$ 分成$\;l$ 个小区间, 对每一个小区间包含的${\tilde x_k}$ 使用中心极限定理[18]构造随机变量${\boldsymbol{S}_i}\tilde x = \frac{1}{T}\sum\nolimits_{k = 1}^T {{{\tilde{\boldsymbol x}}_k}}$ , 得到$\;l$ 个${\boldsymbol{S}_i}\tilde x,\;i \in l$ , 这$l$ 个随机变量是服从标准正态分布的, 当攻击发生时,${S_i}\tilde x$ 会偏离标准整体分布, 从而达到检测攻击的目的.本文所提出的检测方法从系统内部状态误差入手, 不受限于传感器参数的概率分布未知所带来的影响, 不受限于加性噪声和非加性噪声的情况, 可以成功地完成攻击检测, 并从数学角度给出了误检率与阈值之间的关系.
3.1 攻击检测
式(24)定义了随机变量
${\tilde{\boldsymbol x}_k}$ , 本节定义当系统未遭受攻击时, 变量为${\tilde{\boldsymbol x}_k}$ ; 当系统遭受攻击时, 变量为$\tilde{\boldsymbol x}_k^a$ . 下面分别计算这两个随机变量的两个统计特性, 均值和方差.1)当系统正常运行时
$$\begin{split} {\rm{E}}({{\tilde{\boldsymbol x}}_k}) =\;& {\rm{E}}({\boldsymbol{x}_k} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}} - {K_k}({\boldsymbol{y}_k} - {H_k}{{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}))=\\ &{\rm{E}}({\boldsymbol{x}_k} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}} - {K_k}({H_k}{{\hat{\boldsymbol x}}_k} + {v_k} - {H_k}{{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}})) =\\ &{\rm{E}}(({I_k} - {K_k}{H_k})({\boldsymbol{x}_k} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}) - {K_k}{\rm{E}}({v_k}))=\\ &({I_k} - {K_k}{H_k}){\rm{E}}({\boldsymbol{x}_k} - {{\hat{\boldsymbol x}}_{k\left| {k - 1} \right.}}) - {K_k}{\rm{E}}({v_k}) \end{split}$$ (67) 当初值设置为
${\hat{\boldsymbol x}_0} = {\rm{E}}({\boldsymbol{x}_0})$ 时, 式(67)变为$${\rm{E}}({\tilde{\boldsymbol x}_k}) = 0$$ (68) 方差为
$$\begin{aligned} {\hat P}_k|k=\;&{\rm{E}}\left({\tilde{\boldsymbol x}}_{k}\tilde{\boldsymbol x}^{\rm{T}}_k\right)=\\ &{\rm{E}}\left[(I_{k}-K_{k}H_{k})({\boldsymbol{x}}_{k}-\hat{\boldsymbol x}_{k|k-1})-K_{k}{\rm{E}}(v_{k}\right]\times\\ &[(I_{k}-K_{k}H_{k})({\boldsymbol{x}}_{k}-{\hat{\boldsymbol x}}_{k|k-1})-K_{k}{\rm{E}}(v_{k})]^{\rm{T}}=\\ &(I_{k}-K_{k}H_{k})\hat{P}_{k|k-1}(I_{k}-K_{k}H_{k})^{\rm{T}}+K_{k}R_{k}K_{k}^{\rm{T}} \end{aligned}$$ (69) 2)当系统遭受攻击时
$$\begin{split}{\rm{E}}(\tilde{\boldsymbol x}_k^a) = \;&{\rm{E}}(\boldsymbol{x}_k^a - \hat{\boldsymbol x}_{k\left| {k - 1} \right.}^a - {K_k}(\boldsymbol{y}_k^a + {a_k} - {H_k}\hat{\boldsymbol x}_{k\left| {k - 1} \right.}^a))=\\ &{\rm{E}}(\boldsymbol{x}_k^a - \hat{\boldsymbol x}_{k\left| {k - 1} \right.}^a - \\ &{K_k}({H_k}\boldsymbol{x}_k^a + {v_k} + {a_k} - {H_k}\hat{\boldsymbol x}_{k\left| {k - 1} \right.}^a))=\\ &{\rm{ E}}(({I_k} - {K_k}{H_k})(\boldsymbol{x}_k^a - \hat{\boldsymbol x}_{k\left| {k - 1} \right.}^a) - {K_k}{\rm{E}}({v_k})) \end{split}$$ (70) 方差为
$$\begin{split} \hat P_{k\left| k \right.}^a =\;& {\rm{E}}[(\tilde{\boldsymbol x}_k^a){(\tilde{\boldsymbol x}_k^a)^{\rm{T}}}]=\\ & {\rm{E}}\{ [({I_k} - {K_k}{H_k})(\boldsymbol{x}_k^a - \hat{\boldsymbol x}_{k\left| {k - 1} \right.}^a) - {K_k}{\rm{E}}({v_k} + {a_k})] \;\times \\ &{[({I_k} - {K_k}{H_k})(\boldsymbol{x}_k^a - \hat{\boldsymbol x}_{k\left| {k - 1} \right.}^a) - {K_k}{\rm{E}}({v_k} + {a_k})]^{\rm{T}}} =\\ & ({I_k} - {K_k}{H_k})\hat P_{k\left| {k - 1} \right.}^a{({I_k} - {K_k}{H_k})^{\rm{T}}} + {K_k}{R_k}K_k^{\rm{T}} \end{split}$$ (71) 根据第3.1节的描述, 虽然得到了
${\tilde x_k}$ 和$\tilde x_k^a$ 的均值和方差, 但是无法确定随机变量${\tilde x_k}$ 和$\tilde x_k^a$ 服从何种分布. 此时, 对于利用${\tilde x_k}$ 应用中心极限定理的性质, n个独立随机变量之和的分布函数趋于正态分布函数.现构造随机变量
$S\tilde x = \frac{1}{T}\sum\nolimits_{k = 1}^T {{{\tilde{\boldsymbol x}}_k}}$ , 根据中心极限定理, 可知$$S\tilde x \sim {\rm{N}}\left(\frac{1}{T}\sum\limits_{k = 1}^T {{\rm{E}}(\tilde{\boldsymbol x})} ,\frac{1}{T}\sum\limits_{k = 1}^T {{{\hat P}_{k\left| k \right.}}} \right)$$ 同时按照上述形式构造
$S{\tilde x^a} = \frac{1}{T}\sum\nolimits_{k = 1}^T {\tilde{\boldsymbol x}_k^a}$ , 注意此时$S{\tilde x^a}$ 并不服从正态分布, 构造$S{\tilde x^a}$ 是为了保持变量形式上的一致, 记$S{\tilde x^a}$ 服从的分布为$A$ , 得$$S{\tilde x^a} \sim A\left(\frac{1}{T}\sum\limits_{k = 1}^T {{\rm{E}}(\tilde{\boldsymbol x}_k^a)} ,\frac{1}{T}\sum\limits_{k = 1}^T {\hat P_{k\left| k \right.}^a} \right)$$ 为书写简便, 将
$S\tilde x$ 和$S{\tilde x^a}$ 分别记为$$S\tilde x \sim {\rm{N}}(m,\upsilon ),\;\;S{\tilde x^a} \sim A(m + r\upsilon ,q\upsilon )$$ 其中,
$q$ 和$r$ 分别表示系统遭受攻击后均值和方差的变化. 本文使用$m$ 和$\upsilon $ 作为系统正常行为的度量. 接下来, 建立了一对阈值${T_l}$ 和${T_h}$ , 用于检测系统是否遭受攻击,${T_l} = m - k\upsilon $ ,${T_h} = m + k\upsilon $ , 表示允许偏离原数据的最大范围. 当随机变量$S \notin [{T_l},{T_h}]$ 时, 则确定系统遭受攻击[19]. 基于此, 提出一个双边假设检验$$\left\{ {\begin{aligned} {{H_0},\;S \in [{T_l},{T_h}]} \\ {{H_1},\;S \notin [{T_l},{T_h}]} \end{aligned}} \right.\;\;$$ (72) 式中,
${H_0}$ 代表系统没有遭受攻击,${H_1}$ 代表系统遭受攻击.误检率
${P_F}$ 的定义是攻击没有发生却宣布攻击发生的概率, 即$S\tilde x \sim {\rm{N}}(m,\upsilon )$ , 同时,$\left| {S\tilde x - m} \right| > k\upsilon $ .$$\begin{split} {P_F}(k) &= 1 - \int\nolimits_{m - kv}^{m + kv} {\frac{{{\exp (\frac{ - {{(u - m)}^2}}{(2{v^2})})}}}{{v\sqrt {2\pi } }}} {\rm{d}}u =\\ &\;\;\;\;\;1 - \Phi (k) + \Phi ( - k) \\ \end{split} $$ 4. 仿真
本节验证本文设计的自适应平方根无迹卡尔曼滤波算法估计状态及所提出检测算法检测攻击有效性, 并与欧几里得检测方法进行对比. 仿真参数和初始值为
$m = [8.9,{{ }}8.8,{{ }}8.5]$ ,$d = [3.1,{{ }}3.4,{{ }}3.7], $ $u = [6.3, {{ }}1.6,{{ }}8.5]^{\rm{T}}$ ,$B = \left[{0_{3 \times 3}};{\rm{diag}}\{1/m\}\right]$ .为了最小化由噪声引起的误报率, 欧几里得检测法的阈值设计采用
$3\sigma $ 准则($\sigma $ 是噪声信号的标准差), 那么得到先验阈值$f = 3\sigma = 0.3$ .4.1 ASRUKF滤波算法的有效性
对于模型(6), 未遭受攻击时的状态估计, 从图3可以看出, ASRUKF滤波器的滤波性能很好, 滤波精度很高.
4.2 虚假数据注入攻击下的三种检测方法
1) 欧几里得检测法[11]. 由图 4可知, 在迭代步数
$k = 30$ 时刻对系统注入攻击, 在$k = 80$ 时刻成功检测到攻击, 欧几里得检测法可以成功检测出虚假数据注入攻击.2) 巴氏相似性系数检测法[10]. 巴氏系数范围为0到1之间. 越接近0, 证明两序列越相似. 由图5可知, 巴氏相似性系数无法辨别虚假数据注入攻击前后两个残差序列的相似性, 不可以检测出虚假数据注入攻击.
3) 本文所提的检测方法. 首先, 使用分位数–分位数(Quantile-quantile, Q-Q)图来确定前面根据中心极限定理构造的
$S\tilde x$ 分布可以近似为高斯分布. 因为根据中心极限定理, 随机变量的数量越多, 它们的和也越近似服从高斯分布. 然而, 考虑实际情况, 取无限多的独立随机变量是不符合实际的, 所以在近似高斯分布的同时, 要尽量选取较少的独立随机变量. 要利用Q-Q图鉴别样本数据是否近似于高斯分布, 只需观察Q-Q图上的点是否近似地在一条直线附近.本文的检测算法选取了
$T = 30$ 个独立随机变量参与运算, 由图6可知, 数据的分布非常接近直线, 这证明$S\tilde x$ 的分布可以视为高斯分布.由图7可以看出, 本文的检测算法对于隐蔽假数据攻击的有效性, 在迭代步数
$k = 30$ 处注入攻击, 选取30个迭代步骤下的随机变量构造出$S\tilde x$ , 所以在$k = 60$ 处超过了阈值, 实现了攻击检测.图8给出了误检率随
$k$ 变化的曲线.从图8可知, 一个算法的误检率过高, 是因为阈值设置得太低, 导致攻击之外的因素也会使得被检测的量超过阈值. 随着$k$ 值增大, 阈值变高, 误检率自然随之下降.在图8中, 当
$k = {{3}}$ 时, 误检率${P_F} = 0.0027$ , 除此之外, 由图4可知, 欧氏检测法检测到攻击所用时间是$\Delta k = 50$ ; 由图8可知, 检测到攻击所用时间是$\Delta k = 30$ , 检测用时更短.综上, 对比三种检测方法可知, 本文提出的检测方法可以成功检测攻击不受限于传感器参数统计特性未知所带来的影响, 不受限于加性噪声和非加性噪声的情况, 对于实际系统的适用性更强.
5. 结束语
本文研究了智能电网系统中虚假数据注入攻击的检测问题. 针对非线性系统, 噪声统计特性未知的情况, 本文使用自适应平方根无迹卡尔曼滤波算法对系统内部状态和噪声作出估计. 针对传感器参数统计特性未知和非加性噪声的情况, 利用中心极限定理构造出符合正态分布的随机变量, 基于该随机变量提出了一种攻击检测方法, 并从评价指标的角度对算法进行分析, 该算法对于实际系统的适用性更强.
-
[1] 彭大天, 董建敏, 蔡忠闽, 张长青, 彭勤科. 假数据注入攻击下信息物理融合系统的稳定性研究. 自动化学报, 2019, 45(1): 196-205.PENG Da-Tian, DONG Jian-Min, CAI Zhong-Min, ZHANG Chang-Qing, PENG Qin-Ke. On the Stability of Cyber-physical Systems Under False Data Injection Attacks. ACTA AUTOMATICA SINICA, 2019, 45(1): 196-205. [2] 王琦, 邰伟, 汤奕, 倪明. 面向电力信息物理系统的虚假数据注入攻击研究综述. 自动化学报, 2019, 45(1): 72-83.WANG Qi, TAI Wei, TANG Yi, NI Ming. A Review on False Data Injection Attack Toward Cyber-physical Power System. ACTA AUTOMATICA SINICA, 2019, 45(1): 72-83. [3] 童晓阳, 王晓茹. 乌克兰停电事件引起的网络攻击与电网信息安全防范思考. 电力系统自动化, 2016, 40(7): 144-148.Tong Xiao-Yang, Wang Xiao-Ru. Inferrence and Countermeasure Presuppostion of Network Attack in Incdient on Ukrainian Power Grid. Automation of Electric Power Systems, 2016, 40(7): 144-148. [4] Liu Y, Ning P, Reiter M K. False Data Injection Attacks against State Estimation in Electric Power Grids. ACM Transactions on Information and System Security, 2011, 14(1): 1-34. [5] Yang C, Yang W, Shi H. DoS attack in centralised sensor network against state estimation. IET Control Theory & Applications, 2018, 12(9): 1244–1253. [6] Kurt M N, Yılmaz Y, Wang X D. Secure Distributed Dynamic State Estimation in Wide-Area Smart Grids. IEEE Transactions on Information Forensics and Security, 2020, 15: 800-815. doi: 10.1109/TIFS.2019.2928207 [7] Weimer J, Kar S, Johansson K H. Distributed detection and isolation of topology attacks in power networks. In: Proceedings of the 1st ACM International Conference on High Confidence Networked Systems, Beijing, China. 2012. 65−71 [8] Qi J J, Sun K, Wang J H, Liu H. Dynamic State Estimation for Multi-Machine Power System by Unscented Kalman Filter With Enhanced Numerical Stability. IEEE Transactions on smart grid, 2018, 9(2): 1184-1196. doi: 10.1109/TSG.2016.2580584 [9] Mo Y L, Garone E, Casavola A, Sinopoli B. False data injection attacks against state estimation in wireless sensor networks. In: Proceedings of the 49th IEEE Conference on Decision and Control, Atlanta, USA: IEEE, 2010. 5967−5972 [10] Mohammadi A, Plataniotis K N. Noncircular. Attacks on Phasor Measurement Units for State Estimation in Smart Grid. IEEE Journal of Selected Topics in Signal Processing, 2018, 12(4): 777-789. doi: 10.1109/JSTSP.2018.2840517 [11] Manandhar K, Cao X J, Hu F. Detection of Faults and Attacks Including False Data Injection Attack in Smart Grid Using Kalman Filter. IEEE Transactions on Control of Network Systems, 2014, 1(4): 370-379. doi: 10.1109/TCNS.2014.2357531 [12] Yang Q, Yang J, Yu W, An D, Zhang N, Zhao W. On false data-injection attacks against power system state estimation: modeling and countermeasures. IEEE Transactions on Parallel and Distributed Systems, 2014, 25(3): 717–729. doi: 10.1109/TPDS.2013.92 [13] 梁天添, 王茂. 一类广义连续-离散系统量测丢失情况下鲁棒滤波算法. 中国惯性技术学报, 2018, 26(2): 81-86.LIANG Tian-Tian, WANG Mao. Robust algorithm for a class of sampled-data descriptor systems with missing measurements. Journal of Chinese Inertial Technology, 2018, 26(2): 81-86. [14] 赵琳, 王小旭, 孙明, 丁继成, 闫超. 基于极大后验估计和指数加权的自适应UKF滤波算法. 自动化学报, 2010, 36(7): 1008-1019.Zhao Lin, WANG Xiao-Xu, SUN Ming, DING Ji-Cheng, YAN Chao. Adaptive UKF Filtering Algorithm Based on Maximum a Posterior Estimation and Exponential Weighting. ACTA AUTOMATICA SINICA, 2010, 36(7): 1008-1019. [15] 石勇, 韩崇昭. 自适应UKF算法在目标跟踪中的应用. 自动化学报, 2011, 37(6): 755-759.SHI Yong, HAN Chong-Zhao. Adaptive UKF Method with Applications to Target Tracking. ACTA AUTOMATICA SINICA, 2011, 37(6): 755-759. [16] Van der Merwe R, Wan E A. The square-root unscented Kalman filter for state and parameter estimation. In: Proceedings of the 2001 International Conference on Acoustics, Speech, and Signal Processing. New York, USA: IEEE, 2001. 3461−3464 [17] Xiong K, Zhang H Y, Chan C W. Performance Evaluation of UKF-based Nonlinear Filtering. Automatica, 2006, 42(2): 261-270. doi: 10.1016/j.automatica.2005.10.004 [18] Yu W, Griffith D, Ge L Q, Bhattarai S, Golmie N. An Integrated Detection System against False Data Injection Attacks in the Smart Grid. Security and Communication Networks, 2015, 8(2): 91-109. doi: 10.1002/sec.957 [19] Hu J, Wang Z D, Gao H J. Recursive Filtering with Random Parameter Matrices, Multiple Fading Measurements and Correlated Noises. Automatica, 2013, 49(11): 3440-3448. doi: 10.1016/j.automatica.2013.08.021 期刊类型引用(15)
1. 翟千惠，朱萌，俞阳，何玮，康雨萌. 基于XGBoost算法的电力虚假数据注入攻击残差检测. 电子设计工程. 2025(01): 109-112+117 . 
百度学术2. 徐韬，王炜，廖吉星，陈炯. 可拓关联搜索的电力网络错误数据注入检测. 电子设计工程. 2025(02): 154-157+162 . 百度学术3. 黄伟，付红坡，李煜，章卫国. 一种高斯-重尾切换分布鲁棒卡尔曼滤波器. 哈尔滨工业大学学报. 2024(04): 12-23 . 百度学术4. 席磊，田习龙，余涛，程琛. 基于相关特征-多标签级联提升森林的电网虚假数据注入攻击定位检测. 南方电网技术. 2024(05): 39-50+61 . 百度学术5. 王旭，何宇，袁梦薇. 基于Stacking和孤立森林的虚假数据注入攻击防御策略. 智能计算机与应用. 2024(07): 222-226 . 百度学术6. 庞清乐，韩松易，周泰，张峰，焦绪国，王言前. 基于ASRUKF和IMC算法的电力信息物理系统虚假数据注入攻击检测. 智慧电力. 2024(07): 111-118 . 百度学术7. 陈将宏，饶佳黎，李伟亮，胡炀. 基于向量自回归模型的电网虚假数据注入攻击检测. 电力科学与技术学报. 2024(03): 1-9 . 百度学术8. 吴忠强，张伟一. 基于ARO-MKELM的微电网攻击检测. 计量学报. 2024(10): 1444-1452 . 百度学术9. 胡聪，洪德华，张翠翠，王海鑫，薛晓茹，李云路. 一种基于特征映射与深度学习的虚假数据注入检测方法. 现代电力. 2023(01): 125-132 . 百度学术10. 王国庆，杨春雨，马磊，代伟. 基于高斯–广义双曲混合分布的非线性卡尔曼滤波. 自动化学报. 2023(02): 448-460 . 本站查看11. 席磊，何苗，周博奇，李彦营. 基于改进多隐层极限学习机的电网虚假数据注入攻击检测. 自动化学报. 2023(04): 881-890 . 本站查看12. 黄鹏飞，屈剑锋，柴毅，陈小龙，刘切. 强噪声下基于混沌技术的高超声速飞行器传感器故障诊断. 宇航学报. 2023(08): 1203-1212 . 百度学术13. 王新宇，王相杰，张明月，程朋飞，王书征. 基于改进自适应免疫遗传算法的智能电网虚假数据攻击检测方法. 浙江电力. 2023(10): 84-89 . 百度学术14. 杨玉泽，刘文霞，李承泽，刘耕铭，张帅，张艺伟. 面向电力SCADA系统的FDIA检测方法综述. 中国电机工程学报. 2023(22): 8602-8622 . 百度学术15. 夏云舒，王勇，周林，樊汝森. 基于改进生成对抗网络的虚假数据注入攻击检测方法. 电力建设. 2022(03): 58-65 . 百度学术其他类型引用(15)
-
下载:
下载:
计量
- 文章访问数: 1608
- HTML全文浏览量: 511
- PDF下载量: 348
- 被引次数: 30
