2.793

2018影响因子

(CJCR)

  • 中文核心
  • EI
  • 中国科技核心
  • Scopus
  • CSCD
  • 英国科学文摘

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于自适应Kalman滤波的智能电网假数据注入攻击检测

罗小元 潘雪扬 王新宇 关新平

罗小元, 潘雪扬, 王新宇, 关新平. 基于自适应Kalman滤波的智能电网假数据注入攻击检测. 自动化学报, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636
引用本文: 罗小元, 潘雪扬, 王新宇, 关新平. 基于自适应Kalman滤波的智能电网假数据注入攻击检测. 自动化学报, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636
Luo Xiao-Yuan, Pan Xue-Yang, Wang Xin-Yu, Guan Xin-Ping. Detection of False Data Injection Attack in Smart Grid via Adaptive Kalman Filtering . Acta Automatica Sinica, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636
Citation: Luo Xiao-Yuan, Pan Xue-Yang, Wang Xin-Yu, Guan Xin-Ping. Detection of False Data Injection Attack in Smart Grid via Adaptive Kalman Filtering . Acta Automatica Sinica, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636

基于自适应Kalman滤波的智能电网假数据注入攻击检测


DOI: 10.16383/j.aas.c190636
详细信息
    作者简介:

    燕山大学自动化系教授. 2005年获得燕山大学控制科学与工程学科博士学位 主要研究方向为网络控制系统, CPS网络攻击检测等. E-Mail: xyluo@ysu.edu.cn

    燕山大学控制科学与工程专业硕士研究生. 主要研究方向为卡尔曼滤波和智能电网攻击检测. E-Mail: onty123@126.com

    燕山大学控制科学与工程专业博士研究生.主要从事智能电网攻击检测与防御研究. E-Mail: wangxinyuphd@163.com

    上海交通大学自动化系教授. 1999年获得哈尔滨工业大学控制科学与工程学科博士学位.主要研究方向为无线网络系统, CPS网络攻击检测等. E-Mail: xpguan@sjtu.edu.cn

  • 基金项目:  国家自然科学基金(61873228)资助

Detection of False Data Injection Attack in Smart Grid via Adaptive Kalman Filtering

More Information
计量
  • 文章访问数:  3
  • HTML全文浏览量:  4
  • 被引次数: 0
出版历程

基于自适应Kalman滤波的智能电网假数据注入攻击检测

doi: 10.16383/j.aas.c190636
    基金项目:  国家自然科学基金(61873228)资助
    作者简介:

    燕山大学自动化系教授. 2005年获得燕山大学控制科学与工程学科博士学位 主要研究方向为网络控制系统, CPS网络攻击检测等. E-Mail: xyluo@ysu.edu.cn

    燕山大学控制科学与工程专业硕士研究生. 主要研究方向为卡尔曼滤波和智能电网攻击检测. E-Mail: onty123@126.com

    燕山大学控制科学与工程专业博士研究生.主要从事智能电网攻击检测与防御研究. E-Mail: wangxinyuphd@163.com

    上海交通大学自动化系教授. 1999年获得哈尔滨工业大学控制科学与工程学科博士学位.主要研究方向为无线网络系统, CPS网络攻击检测等. E-Mail: xpguan@sjtu.edu.cn

摘要: 本文研究了一种针对智能电网中假数据注入攻击的有效检测方法. 假数据注入攻击可以保持攻击前后残差基本不变, 绕过传统的不良数据检测技术. 首先基于电网模型, 分析了假数据注入攻击的攻击特性, 针对噪声统计特性未知且无迹Kalman滤波不稳定的现象, 提出了自适应平方根无迹Kalman滤波改进算法. 基于状态估计值, 结合中心极限定理提出检测算法, 并与欧几里得检测方法, 巴氏系数检测方法作比较. 最后, 仿真表明本文所提检测算法的优越性.

English Abstract

罗小元, 潘雪扬, 王新宇, 关新平. 基于自适应Kalman滤波的智能电网假数据注入攻击检测. 自动化学报, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636
引用本文: 罗小元, 潘雪扬, 王新宇, 关新平. 基于自适应Kalman滤波的智能电网假数据注入攻击检测. 自动化学报, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636
Luo Xiao-Yuan, Pan Xue-Yang, Wang Xin-Yu, Guan Xin-Ping. Detection of False Data Injection Attack in Smart Grid via Adaptive Kalman Filtering . Acta Automatica Sinica, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636
Citation: Luo Xiao-Yuan, Pan Xue-Yang, Wang Xin-Yu, Guan Xin-Ping. Detection of False Data Injection Attack in Smart Grid via Adaptive Kalman Filtering . Acta Automatica Sinica, 2020, 41(x): 1−12. doi: 10.16383/j.aas.c190636
  • 智能电网是一种新型的电网, 它采用先进的通信网络技术和控制技术来支持更高效的能源安全传输和分配. 然而, 由于智能电网系统的复杂性和开放性, 智能电网中进行数据交换的信息网络成为易受到恶意攻击的对象[1-2]. 例如, 2016年, 黑客攻击乌克兰国家电力部门致使国内发生了一次大规模的停电事件[3], 造成严重经济损失. 因此, 智能电网的攻击检测研究具有重要意义.

    隐蔽假数据攻击是目前恶意攻击的典型代表, 攻击者对传感器节点注入精心设计的错误数据, 接收错误数据的控制中心, 继而做出错误决策破坏系统的稳定性[4]. 在文献[5]中, 拒绝服务攻击旨在中断电力网络通信信道的可用性. 文献[6]构建了一种对智能电网中以完整性和可用性为目标的攻击分类方法. 文献[7]解决了在电力网络中攻击检测和拓扑隔离的问题. 文献[8]提出了一种算法来识别要操作的智能电表的最优数目, 从而找到最优攻击策略, 其目的是干扰电网系统的状态估计, 影响其稳定性. 文献[9]设计了具有隐蔽特性的虚假数据攻击, 它可以使攻击前后残差基本不变, 因此, 基于卡方检测器的检测技术是无效的. 近年来, 隐蔽性攻击检测成为了研究热点之一.

    针对智能电网遭受虚假数据注入攻击的检测问题, 近年来有了很多成果. 文献[10]中, 以电机的电压模型为研究对象, 使用卡尔曼滤波技术得到该系统的残差序列, 通过计算攻击前后残差序列之间的Bhattacharyya距离来判断系统中是否存在攻击. 该文献的不足之处有两点: 其一, Bhattacharyya距离无法辨别虚假数据注入攻击前后两个残差序列的相似性, 因为虚假数据注入攻击前后残差保持不变; 其二, 该电压模型是线性的, 对于实际中存在的大多数的非线性系统, 该检测算法是失效的. 文献[11]使用卡尔曼滤波技术得到系统的状态, 从系统状态角度考虑设计检测算法, 提出欧几里德检测方法, 该方法可以检测隐蔽假数据注入攻击. 该文献的不足之处是, 系统方程对于噪声是线性的, 即该噪声是加性噪声, 当系统方程对于噪声是非线性的时候, 噪声经过非线性变换, 不再服从高斯分布, 也就无法设计阈值, 无法检测攻击. 文献[12]中, 攻击者可能会缓慢改变多个传感器的测量值, 因此上述统计异常检测不会检测到个别受损的测量值, 所以提出检测思想, 这些测量值组合起来会导致状态变量远离其真值, 然而, 文中作出假设, 控制中心收集到的测量值都是服从高斯分布的, 并基于此性质设计了双边假设检验检测攻击. 该文献的不足之处是, 电网系统中的参数仅会在特定情况下服从高斯分布, 而这种情况并不常见, 且建模时没有考虑噪声, 所以该文中的检测算法局限性较大.

    因此, 根据虚假数据注入攻击的特性, 考虑系统非线性和噪声统计特性未知情况, 本文提出一种基于自适应平方根无迹卡尔曼滤波器的智能电网隐蔽假数据攻击检测方法. 该算法可以有效地, 稳定地应用于非线性系统当中对系统状态作出估计, 依据系统状态设计检测算法检测攻击, 并从检测指标的角度与现有算法作对比. 最后进行仿真实验, 实验结果证明所提出的基于该算法的攻击检测方法可以准确给出状态估计值从而检测出隐蔽攻击.

    • 我们考虑一个由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)

      图  1  3总线电网模型

      Figure 1.  3-bus grid model

      $${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{gathered} {{\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{gathered} $$ (2)

      每一条总线都安装有传感器, 将总线的数据传输回控制中心, 系统输出方程如下

      $${{y}_{{i}}} = \left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right]{{x}_{{i}}} + {\rho }(t)$$ (3)

      ${{y}_{{i}}}$表示总线${i}$的输出向量, ${\rho }$是测量噪声.

      针对上述系统, 本文采用欧拉离散化方法[13], 即对于一个非线性系统

      $$\dot{ x} = f({x}) + B\tau + \varpi $$

      离散化后的形式,

      $${{x}_{{k}}} = {{x}_{{k} - 1}} + \tau f({{x}_{{k} - 1}}) + \tau B{u} + \tau \varpi + {\rm{O}}({\tau ^2})$$

      其中

      $${\rm{O}}({\tau ^2}) = \dfrac{{{\tau ^2}}}{2}\dot f({x}) + \dfrac{{{\tau ^3}}}{3}{f^{(2)}}({x}) + \cdots + \dfrac{{{\tau ^{\rm{k}}}}}{{k!}}{f^{({k} - 1)}}({x})$$

      采样间隔值$\tau $取值充分小以保证高阶项${\rm{O}}({\tau ^2})$忽略不计, 本文中$\tau = 0.1s$.

      对(2)使用上述欧拉离散化算法对该系统离散化, 得到下式,

      $$\begin{array}{l} {{x}_{{k}}} = {f_{{o}}}({{x}_{{k} - 1}}) + {B_o}{u} + {{\eta }_{k - 1}}\\ {{y}_{{k}}} = H{{x}_{{k}}} + {{\rho }_k} \end{array}$$ (4)

      其中, ${f_o}( \cdot )$表示关于$x$的非线性函数, 将${{x}_{{k} - 1}} + \tau {{f}}({{x}_{{k} - 1}})$合并在一起, ${B_o} = \tau B$, $\eta = \tau \varpi $.

      假设$\eta $$\rho $是互不相关的非零均值高斯白噪声, 统计特性如下,

      $$\left\{ {\begin{array}{*{20}{l}} {E({{{{\eta }}}_k}) = {{{q}}},{{ }}{\rm{cov}} ({{{{\eta }}}_k}{{{\eta }}}_j^T) = {Q_k}{\delta _{kj}}} \\ {E({{{{\rho }}}_k}) = {{{r}}},{{ }}{\rm{cov}} ({{{{\rho }}}_k}{{{\rho }}}_j^T) = {R_k}{\delta _{kj}}} \\ {E({\eta _k}\rho _j^T) = 0} \end{array}} \right.$$ (5)

      式中, Q是非负定矩阵, R是正定矩阵, ${\delta _{kj}}$$Krone cher - \delta$函数,

      $${\delta _{kj}} = \left\{ {\begin{array}{*{20}{c}} {0{{ }}if{{ }}k \ne j} \\ {1{{ }}if{{ }}k = j} \end{array}} \right.$$

      $\eta $$\rho $改写为${\eta _{{k}}} = {q} + {\mu _{{k}}}$, ${\rho _{{k}}} = {r} + {{v}_{{k}}}$, 有

      $$\left\{ {\begin{array}{*{20}{l}} {{{x}_{{k}}} = {f_o}({{x}_{{k} - 1}}) + {B_o}{u} + {q} + {{\rho }_{k - 1}}}\\ {{{y}_{{k}}} = H{{x}_{{k}}} + {r} + {{v}_{k - 1}}} \end{array}} \right.$$ (6)

      ${\mu _k}$${v_k}$是互不相关的零均值高斯白噪声.${{x}_{{k}}} = {[{\delta _1}{\delta _2}{\delta _3}{\omega _1}{\omega _2}{\omega _3}]^{\rm{T}}}$, ${{y}_{{k}}} \in {\bf{R}^6}$.

    • 假设一个恶意的第三方想要破坏第1节中描述的系统的完整性. 本文主要考虑虚假数据注入攻击.

      假设攻击者具有系统知识, 知道系统矩阵, 控制矩阵, 测量矩阵, 可以控制一系列系统中的传感器数据.

      考虑假数据注入攻击模型描述如下

      $${{z}_{{a}}}(k) = \left\{ {\begin{array}{*{20}{l}} {H{{x}_{{k}}} + \xi (k)k\not \subset T}\\ {H{{x}_{{a,k}}} + \xi (k) + {{y}_{{a}}}(k)k \subset T} \end{array}} \right.$$ (7)

      式中, $\varGamma = {{\rm{diag}}}({\gamma _1},...,{\gamma _n})$代表传感器选择矩阵, T是攻击的时间范围, 当${\gamma _i} = 1$时, 代表第$i$个传感器遭受攻击, 否则${\gamma _i} = 0$. 且${{y}_{{a}}}(k)$是攻击者精心构建的攻击向量序列, 攻击框图如图2所示.

      图  2  系统遭受攻击框图

      Figure 2.  block diagram of System under attack

      在隐蔽假数据注入攻击情况下, 在所考虑的攻击模型中, 攻击者构造攻击序列注入传感器, 在该序列下, 状态差值最终将发散到$\infty $, 而不会触发卡方检测器. 该攻击序列满足以下性质[9]:

      $$\begin{split} &\mathop {\lim }\limits_{k \to T} \left\| {{{x}_{{a}}}(k) - {x}(k)} \right\| \geqslant \alpha (\alpha > 0),\\ &\mathop {\lim }\limits_{k \to T} \left\| {{{r}_{{a}}}(k)} \right\| \leqslant \varepsilon (0 < \varepsilon < 1). \end{split}$$ (8)

      其中, ${{x}_{{a}}}({{k}})$${{r}_{{a}}}({{k}})$是遭受攻击时系统的状态变量和残差. 对攻击者克服检测机制的隐蔽性进行分析[4]. 由输出方程可得

      $$\begin{split} {{r}_a} =& \left\| {({y} + {a}) - H(\hat{ x} + {c})} \right\|\\ =& \left\| {({y} - H\hat{ x}) + ({a} - H{c})} \right\|\\ = & \left\| {{y} - H\hat{ x}} \right\| + \left\| {{a} - H{c}} \right\| = {r} + {\tau _a} \end{split}$$ (9)

      式中${a}$指攻击信号, $\hat{ x}$指状态估计值, ${c}$指状态变化量. 由式(9)知, 当${a} = H{c}$时, ${{r}_{{a}}} = {r}$, 从输出残差角度考虑, 实现了隐蔽性.

      由上述分析可知, 从输出残差角度无法实现攻击检测. 所以, 需要借助滤波器从系统状态入手, 利用状态差值检测攻击信号. 本文考虑通过基于自适应平方根无迹卡尔曼滤波器的方法, 从系统状态入手来检测此种攻击.

    • 本文提出将噪声估计环节[14]加入到无迹卡尔曼滤波算法中, 实现对于噪声统计特性的在线估计, 同时改变标准UKF中状态误差协方差矩阵的迭代方式来保证滤波器的稳定性.

    • 本小节给出平方根无迹变换的实现方式.

      步骤 1. 利用${{k}} - 1$时刻的估计状态及状态误差协方差矩阵来计算sigma采样点$\xi {}_i{{ }}$并给出其权值${W_{m,i}}$, ${W_{c,i}}$.

      $$\left\{ {\begin{array}{*{20}{l}} {{\xi _{{i},k - 1}} = {{\hat{ x}}_{k - 1}},i = 1}\\ {{\xi _{{i},k - 1}} = {{\hat{ x}}_{k - 1}} + \gamma {{({S_{k - 1}})}_i},i = 2, \ldots ,n + 1}\\ {{\xi _{{i},k - 1}} = {{\hat{ x}}_{k - 1}} - \gamma {{({S_{k - 1}})}_i},i = n + 2, \ldots ,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} )^T}$, 故${({{S}_{k - 1}})_i}$等于$(\sqrt {(n + \lambda )P} )_i^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}} \leqslant \alpha \leqslant 1$, $\beta = 2$, $n + \kappa = 3$[8].

      步骤 2.计算sigma点通过非线性函数的传播结果

      $$\begin{split} &{\tau _i} = f({\xi _i})i = 1, \ldots ,2n + 1\\ &{z} = \sum\limits_{i = 1}^{2n + 1} {{W_{m,i}}{{\tau }_i}} \\ &{P_x} = \sum\limits_{i = 1}^{2n + 1} {{W_{c,i}}({{\tau }_i} - {z}){{({{\tau }_i} - {z})}^{\rm{T}}}} \end{split}$$ (10)

      式(10)是无迹变换中的计算误差协方差矩阵的公式, 平方根无迹变换把(11)换为

      $$\begin{split}{S_z} =& {qr}(\sqrt {{{W}_{c,1}}} ({\tau _1} - {z}) \ldots \\ &\sqrt {{{W}_{c,2{n} + 1}}} ({\tau _{2n + 1}} - {z}))\end{split}$$ (11)
    • (1)假设噪声是非零均值的高斯噪声, 本文使用的是Sage-Husa噪声估计器, 通过极大后验估计原理, 获得次优噪声估计值, 噪声估计部分参考了文献[14].

      对于系统(6), 系统方程是非线性的, 输出方程是线性的, 所以噪声估计器如下,

      $${\hat{ q}_{k + 1}} = \frac{1}{k}\left[(k - 1){\hat{ q}_k} + {\hat{ 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}}{{\varepsilon }_{k + 1}}{{\varepsilon }_{k + 1}}K_{k + 1}^{\rm{T}} \!+\! {P_{k + 1}}-\\ & \sum\limits_{i = 1}^{2n + 1} {W_i^c} ({{\tau }_{i,k + 1\left| k \right.}} - {{\hat{ x}}_{k + 1\left| k \right.}})\\ &{({{\tau }_{i,k + 1\left| k \right.}} - {{\hat{ x}}_{k + 1\left| k \right.}})^{\rm{T}}}\big] \end{split}$$
      $$\begin{split} &{{\hat{ r}}_{k + 1}} = \dfrac{1}{k}\left[(k - 1){{\hat{ r}}_k} + {{y}_{k + 1}} - H{{\hat{ x}}_{k + 1\left| k \right.}}\right]\\ &{{\hat{ R}}_{k + 1}} = \dfrac{1}{k}\left[(k - 1){{\hat{ R}}_k} + {{\varepsilon }_{k + 1}}{\varepsilon }_{k + 1}^{\rm{T}} - H{P_{k + 1\left| k \right.}}H_k^{\rm{T}}\right] \end{split}$$

      ${{\varepsilon }_k} = {{y}_k} - {\hat{ y}_{k\left| {k - 1} \right.}}$代表残差, $K$代表滤波器增益矩阵.

    • 本节用到了MATLAB中的3个函数, $qr( \cdot )$表示QR分解, $A \in {\mathbb{R}^{m \times l}}$, $[Q,R] = $ $qr(A)$, 该函数生成$m \times l$上三角形矩阵R和$m \times m$酉矩阵Q, 这样$A = Q*R$.

      $R{{1}} = cholupdate(R,x)$, 返回$A + xx'$的上三角Cholesky因子, $x$是具有合适长度的一个列向量, 其中$R = chol(A)$A的原始Cholesky分解因子.

      $diag( \cdot )$函数, $D = diag(v)$, 返回以向量v的元素为主对角线的对角矩阵.$x = diag(A)$, 返回A的主对角线元素的列向量. 对于式(22)中最外侧的$diag( \cdot )$, 返回一个对角矩阵, 主对角线元素是每一次迭代得到的噪声矩阵主对角线元素的平方根的实时估计值.

      自适应平方根无迹卡尔曼滤波算法的具体步骤如下

      (1) 初始化(k=1):

      $$\left\{ {\begin{array}{*{20}{l}} {{{\hat{ x}}_1} = E[{{x}_1}]}\\ {{S_1} = chol\left\{ E\left[({{x}_1} - {{\hat{ x}}_1}){{({{x}_1} - {{\hat{ x}}_1})}^{\rm{T}}}\right]\right\} }\\ {\sqrt {{{\hat Q}_1}} = {S_1}}\\ {{{\hat{ q}}_1} = 0} \end{array}} \right.$$ (12)

      (2) 对于k=2,3... 进行迭代

      时间更新

      计算sigma点, 构造矩阵

      $${{\xi }_{k - 1}} = \left[{\hat{ x}_{k - 1}}{\hat{ x}_{k - 1}} + \gamma {({S_{k - 1}})_i}{\hat{ x}_{k - 1}} - \gamma {({S_{k - 1}})_i}\right]$$ (13)
      $$\left\{ {\begin{array}{*{20}{l}} {{\tau _{i,k\left| {k - 1} \right.}} = f({{\xi }_{i,k - 1}}) + {q}}\\ {{x} = \sum\limits_{i = 1}^{2n + 1} {{W_{m,i}}{\tau _{i,k\left| {k - 1} \right.}} + \hat{ q}} }\\ {{S_{k\left| {k - 1} \right.}}\! =\! qr\left\{ \!\left[\!\sqrt {{W_{c,i}}} ({\tau _{2:2n + 1}} - {{\hat{ x}}_{k\left| {k - 1} \right.}})\sqrt {{{\hat Q}_{k - 1}}} \right]\right\} }\\ {S = cholupdate}\\ {\left\{ \left[{S_{k\left| {k - 1} \right.}},\sqrt {{W_{c,i}}} ({\tau _{2:2n + 1}} - {{\hat{ x}}_{k\left| {k - 1} \right.}}), \pm \right]\right\} } \end{array}} \right.$$ (14)

      测量更新

      $${\tilde{ y}_k} = {{y}_k} - H{\hat{ x}_{k\left| {k - 1} \right.}}$$ (15)

      计算混合误差协方差矩阵平方根

      $${P_{xy}} = S_{k\left| {k - 1} \right.}^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]^T}\bigg\} $$ (17)

      计算滤波增益

      $${K_k} = \left({P_{xy,k}}/S_{y,k}^T\right)/{S_{y,k}}$$ (18)

      更新状态值

      $${\hat{ x}_k} = {\hat{ x}_{k\left| {k - 1} \right.}} + {K_k}{\tilde{ y}_k}$$ (19)

      更新状态误差协方差矩阵平方根

      $$\left\{ {\begin{array}{*{20}{l}} {U = {K_k}S_{y,k}^T} \\ {{S_k} = cholupdate\left\{ {S_{k\left| {k - 1} \right.}}{{ }}U{{ }} - 1\right\} } \end{array}} \right.$$ (20)

      估计过程噪声均值

      $${\hat{ q}_{k + 1}} = \dfrac{1}{k}\left[(k - 1){\hat{ q}_k} + {\hat{ x}_{k + 1}} - \sum\limits_{i = 1}^{2n + 1} {W_i^m} {f_o}({{\xi }_{i,k}})\right]$$ (21)

      估计过程噪声协方差矩阵平方根

      $$\left\{ {\begin{array}{*{20}{l}} {\sqrt {Q1} = cholupdate\left\{ \sqrt {{{\hat Q}_{k - 1}}} ,\left| {{{\hat{ x}}_k} - {{\hat{ 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}} = diag\left\{ \sqrt {diag\left(\sqrt {Q2} {{\sqrt {Q2} }^T}\right)} \right\} }\\[-10pt] \end{array}} \right.$$ (22)

      算法结束.

      注1: 1.Sage-Husa噪声估计器不能同时处理过程噪声和测量噪声都未知的情况, 会造成滤波发散, 上述算法假设测量噪声的均值和方差是已知的[15].

      注2: 式(14)最后一行中, $ \pm $表示当${W_{c,1}} > 0$时, 取正, 反之, 取负, 为保证根号下不出现负数[16].

    • 本小节给出自适应平方根无迹卡尔曼滤波算法的估计误差保持随机有界性的充分条件并给出证明.

      对于非线性电网系统(6)

      $$\begin{split} &{{x}_k} = {f_o}({{x}_{k - 1}}) + {B_o}{u} + {q} + {{\mu }_{k - 1}}\\ &{{y}_k} = H{{x}_{k - 1}} + {r} + {{v}_k} \end{split}$$ (23)

      定义

      $${\tilde{ x}_k} = {{x}_k} - {\hat{ x}_k}$$ (24)
      $${\tilde{ x}_{k\left| {k - 1} \right.}} = {{x}_k} - {\hat{ x}_{k\left| {k - 1} \right.}}$$ (25)

      ${\hat{ x}_{k - 1}}$处对${{{{x}}}_k}$泰勒级数展开

      $$\begin{split} &{{x}_k} = {f_o}({{\hat{ x}}_{k - 1}}) + {B_o}{u} + {q} + \\ &\nabla {f_o}({{\hat{ x}}_{k - 1}}){{\tilde{ x}}_{k - 1}}\! +\! \dfrac{1}{2}{\nabla ^2}{f_o}({{\hat{ x}}_{k - 1}})\tilde{ x}_{k - 1}^2 \!+\! \ldots \!+\! \omega \end{split}$$ (26)

      其中

      $${\nabla ^i}{f_{{o}}}(\hat{ x}){\tilde{ x}^i} = {\left(\sum\limits_{j = 1}^n {{{\tilde{ x}}_j}} \frac{\partial }{{\partial {{x}_j}}}\right)^i}{\left. {{f_{{o}}}({x})} \right|_{{x} = {{\hat{ x}}_{k - 1}}}}$$

      ${\hat x_{k - 1}}$处对${\hat x_{k\left| {k - 1} \right.}}$泰勒级数展开[19]

      $$\begin{split} {{\hat{ x}}_{k\left| {k - 1} \right.}} = &\left(1 - \frac{n}{{n + \lambda }}\right) \times (f({{\hat{ x}}_{k - 1}}) + B{u}) + {q} + \\ &\frac{1}{{2(n + \lambda )}} \times \sum\limits_{i = 2}^{n + 1} {(f\left[({{\hat{ x}}_{k - 1}}) + \gamma {{({S_{k - 1}})}_i}\right] + } \\ &B{u} + {q}) + \frac{1}{{2(n + \lambda )}} \times \sum\limits_{i = 2}^{n + 1} {\bigg(f\bigg[({{\hat{ x}}_{k - 1}}) - } \\ &\gamma {({S_{k - 1}})_i}\bigg] + B{u} + {q}\bigg)\\ =& f({{\hat{ x}}_{k - 1}}) + B{u} + {q} + {\nabla ^2}f({{\hat{ x}}_{k - 1}}){P_{k - 1}}\\[-10pt] \end{split}$$ (27)

      将式(26)和式(27)代入式(25)

      $${\tilde{ x}_{k\left| {k - 1} \right.}} = {\beta _k}{F_k}{\tilde{ x}_{k - 1}} + {{\mu }_k}$$ (28)

      式中

      $${F_k} = {\left. {\dfrac{{\partial {f_o}({{x}_{k - 1}})}}{{\partial {{x}_{k - 1}}}}} \right|_{{x} = {{\hat{ x}}_{k - 1}}}}$$ (29)

      ${\beta _k} = diag\left\{ {{\beta _{1k}}, \cdots ,{\beta _{nk}}} \right\}$为未知的矩阵, 其作用是弥补建模一阶线性化引起的模型误差.

      对于(24)中的输出方程, 定义

      $${\tilde{ y}_{k\left| {k - 1} \right.}} = {{y}_k} - {\hat{ y}_{k\left| {k - 1} \right.}}$$ (30)

      式(30)可写成

      $$\begin{split} {{\tilde{ y}}_{k\left| {k - 1} \right.}} =& H{{x}_k} + {r} + {{v}_k} - (H{{\hat{ x}}_{k\left| {k - 1} \right.}} + {r})\\ = &H{{\tilde{ x}}_{k\left| {k - 1} \right.}} + {r} + {{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({{\tau }_{2:2n + 1,k\left| {k - 1} \right.}} - {\hat{ x}_{k\left| {k - 1} \right.}}\right)\sqrt {{{\hat Q}_{k - 1}}} \right]^{\rm{T}}}$$

      使用$qr( \cdot )$函数有$A = Q*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({{\tau }_{i,k\left| {k - 1} \right.}} - {{\hat{ x}}_{k\left| {k - 1} \right.}}\right)}^{\rm{T}}}\\ &\left({{\tau }_{i,k\left| {k - 1} \right.}} - {{\hat{ 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({{\tau }_{1,k\left| {k - 1} \right.}} \!-\! {{\hat{ x}}_{k\left| {k - 1} \right.}}\right) \!\times\! \\ &\sqrt {{W_{c,1}}} {\left({{\tau }_{1,k\left| {k - 1} \right.}} - {{\hat{ 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.}} = &E\left[{{\tilde{ x}}_{k\left| {k - 1} \right.}}\tilde{ x}_{k\left| {k - 1} \right.}^{\rm{T}}\right]\\ =& E\left[({\beta _k}{F_k}{{\hat{ x}}_{k - 1}} + {{\mu }_k}){({\beta _k}{F_k}{{\hat{ x}}_{k - 1}} + {{\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.}}$$E[({\beta _k}{F_k}{\tilde{ x}_{k - 1}}){({\beta _k}{F_k}{\tilde{ x}_{k - 1}})^{\rm{T}}}]$${\beta _k}{F_k} {\hat P_{k - 1}} F_k^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^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^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^T} + {R_k}$$ (35)
      $${\hat P_{xy}} = {\hat P_{k\left| {k - 1} \right.}}{H^T}$$ (36)
      $${\hat P_k} = {\hat P_{k\left| {k - 1} \right.}} - {\hat P_{xy}}\hat P_{yy}^{ - 1}\hat P_{xy}^T$$ (37)

      引理 1[17] 对于一个随机变量${\xi _k}$及其函数$V({\xi _k})$, 实数${\upsilon _{\min }},{\upsilon _{\min }},\mu > 0,0 < \lambda \leqslant 1$对于任意的k, 若满足以下两个不等式性质,

      $$\begin{split} &{\upsilon _{\min }}{\left\| {{{\xi }_k}} \right\|^2} \leqslant V({\xi _k}) \leqslant {\upsilon _{\max }}{\left\| {{{\xi }_k}} \right\|^2},\\ &E\left[V({{\xi }_k})\left| {{{\xi }_{k - 1}}} \right.\right] - V\left({{\xi }_{k - 1}}\right) \leqslant \mu - \lambda V\left({{\xi }_{k - 1}}\right) \end{split}$$ (38)

      则该随机变量均方有界, 即

      $$\begin{split} E[{\left\| {{\xi _k}} \right\|^2}] \leqslant& \dfrac{{{\upsilon _{\max }}}}{{{\upsilon _{\min }}}}E[{\left\| {{\xi _0}} \right\|^2}]{(1 - \lambda )^k} + \\ &\dfrac{\mu }{{{\upsilon _{\min }}}}\sum\limits_{i = 1}^{k - 1} {{{(1 - \lambda )}^i}.} \\ \end{split} $$ (39)

      引理 2[17] 假设矩阵$A \in {R^{n \times m}}$, $B \in {R^{m \times n}}$, $C \in {R^{n \times n}}$, 如果$A > 0,C > 0$, 则有

      $${A^{ - 1}} > B{({B^T}AB + C)^{ - 1}}{B^T}$$ (40)

      引理 3[17] 假设矩阵$A \in {R^{n \times n}}$, $C \in {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 \leqslant &{F_k}F_k^{\rm{T}} \leqslant f_{\max }^2I \\ h_{\min }^2I \leqslant& {H_k}H_k^{\rm{T}} \leqslant h_{\max }^2I \\ \beta _{\min }^2I \leqslant &{\beta _k}\beta _k^{\rm{T}} \leqslant \beta _{\max }^2I \\ \end{split} $$ (42)

      (2) 存在实数${\Xi _{\max }},{\Xi _{\min }},{r_{\min }},{p_{\max }},{p_{\min }}$, 使得下式成立:

      $$\begin{split} & {p_{\min }}I \leqslant {P_k} \leqslant {p_{m{{ax}}}}I \\ & {\Xi _{\min }}I < {\Xi _k} \leqslant {\Xi _{\max }}I \\ &{r_{\min }}I \leqslant {R_k} \\ \end{split} $$ (43)

      那么自适应平方根无迹卡尔曼滤波算法的状态估计误差将是均方有界的, 即自适应平方根无迹卡尔曼滤波器稳定.

      证明 : 选择函数

      $${V_k}({\tilde x_k}) = \tilde x_k^T\hat P_k^{ - 1}{\tilde x_k}$$ (44)

      由式(43)可知

      $$\dfrac{1}{{{p_{\max }}}}{\left\| {{{\tilde x}_k}} \right\|^2} \leqslant V({\tilde x_k}) \leqslant \dfrac{1}{{{p_{\min }}}}{\left\| {{{\tilde x}_k}} \right\|^2}$$

      满足了引理1中第一个性质, 为满足引理1, 还需要得到式(38)中第二个性质$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}^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^T{({H_k}{\hat P_{k\left| {k - 1} \right.}}H_k^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{ x}}_k}) =& {\left({{\tilde{ x}}_{k\left| {k - 1} \right.}} - {K_k}{{\tilde{ y}}_k}\right)^T}\hat P_k^{ - 1}\left({{\tilde{ x}}_{k\left| {k - 1} \right.}} - {K_k}{{\tilde{ y}}_k}\right)\\ =& {{\tilde{ x}}_{k\left| {k - 1} \right.}}\hat P_k^{ - 1}{{\tilde{ x}}_{k\left| {k - 1} \right.}} - {\left[{H_k}{{\tilde{ x}}_{k\left| {k - 1} \right.}} + {v_k}\right]^T} \times \\ &K_k^T\hat P_k^{ - 1}{{\tilde{ x}}_{k\left| {k - 1} \right.}} \!-\! \tilde x_{k\left| {k - 1} \right.}^T\hat P_k^{ - 1}{K_k}\big[{H_k}{{\tilde x}_{k\left| {k - 1} \right.}} \!+\! \\ &{v_k}\big] + {\left[{H_k}{{\tilde{ x}}_{k\left| {k - 1} \right.}} + {v_k}\right]^T}K_k^T \times \\ &\hat P_k^{ - 1}{K_k}\left[{H_k}{{\tilde{ 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^TR_k^{ - 1} = {\hat P_k}H_k^TR_k^{ - 1}$$ (49)

      式(45)用矩阵求逆公式得

      $$\hat P_k^{ - 1} = \hat P_{k\left| {k - 1} \right.}^{ - 1} + H_k^TR_k^{ - 1}{H_k}$$ (50)

      将式(49)-(50), 式(28)代入式(48)中, 式(48)的条件期望

      $$\begin{split} E[{V_k}({{\tilde x}_k})\left| {{{\tilde x}_{k\! -\! 1}}} \right.] \!=\!& E{\bigg\{} {({A_{k - 1}}{{\tilde x}_{k - 1}} + {\mu _k})^T}\\ &\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]^T}\\ &\left(\Sigma _k^{ - 1} - \Sigma _k^{ - 1}{B_k}{{\hat P}_k}B_k^T\Sigma _k^{ - 1}\right) \times \\ &\left[{B_k}({A_{k - 1}}{{\tilde x}_{k - 1}} + {\mu _k})\right] + \\ &v_k^T\Sigma _k^{ - 1}{B_k}{{\hat P}_k}B_k^T\Sigma _k^{ - 1}{v_k}\left| {{{\tilde x}_{k - 1}}} \right.\!{\bigg\}} \end{split}$$ (51)

      将式(34)代入式(51)中, 并由引理3得

      $$\begin{split} &{E[{V_k}({{\tilde{ x}}_k})\left| {{{\tilde{ x}}_{k - 1}}} \right.]} \leqslant\\ &{ E\left\{ {{{({\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})}^T}} \right.{{\left( {{\beta _k}{F_k}\hat P_{k{\rm{ - }}1}^{}F_k^T{\beta _k}} \right)}^{ - 1}}}\\ &{({\beta _k}{F_k}{{\tilde{ x}}_{k - 1}}) + \omega _k^T{{({\beta _k}{F_k}\hat P_{k{\rm{ - }}1}^{}F_k^T{\beta _k} + {{\hat Q}_k})}^{ - 1}} \times }\\ &{{\omega _k} - {{({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})}^T} \times \left[ {{H_k}} \right.\left( {{\beta _k}} \right.{F_k}{{\hat P}_{k - 1}}F_k^T \times }\\ &{{\beta _k} + \left. {{{\hat Q}_k}} \right)H_k^T + {{\left. {{R_k}} \right]}^{ - 1}}({H_k}{\beta _k}{F_k} \times {{\tilde{ x}}_{k - 1}}) - }\\ &{{{({H_k}{\omega _k})}^T}\left[ {{H_k}} \right.\left({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k}\right) \times H_k^T + }\\ &{{{\left. {{R_k}} \right]}^{ - 1}}({H_k}{\omega _k}) + v_k^TR_k^{ - 1}{H_k}{{\hat P}_k}H_k^TR_k^{ - 1}{v_k}\left. {\left| {{{\tilde{ x}}_{k - 1}}} \right.} \right\}} \end{split}$$ (52)

      观察式(52)的第一项

      $$\begin{split} &E\bigg\{ {({\beta _k}{F_k}{{\tilde x}_{k - 1}})^T}{\left( {{\beta _k}{F_k}\hat P_{k{\rm{ - }}1}^{}F_k^T{\beta _k}} \right)^{ - 1}} \times \\ &({\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})\left| {{{\tilde{ x}}_{k - 1}}} \right.\bigg\} \\ &= \tilde{ x}_{k - 1}^T{{\hat P}_{k - 1}}{{\tilde{ x}}_{k - 1}} = {V_{k - 1}}({{\tilde{ x}}_{k - 1}}) \end{split}$$ (53)

      两边同时减去式(53)

      $$\begin{split} &{E[{V_k}({{\tilde{ x}}_k})\left| {{{\tilde{ x}}_{k - 1}}} \right.] - {V_{k - 1}}({{\tilde x}_{k - 1}})}\leqslant\\ &{ E\left\{ {\omega _k^T} \right.{{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})}^{ - 1}} - }\\ &{{{({H_k}{\omega _k})}^T} \times \left[ {{H_k}} \right.\left( {{\beta _k}} \right.{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + }\\ &{\left. {{{\hat Q}_k}} \right)H_k^T + {{\left. {{R_k}} \right]}^{ - 1}}({H_k}{\omega _k}) + v_k^TR_k^{ - 1}{H_k} \times }\\ &{{{\hat P}_k}H_k^TR_k^{ - 1}{v_k}\left. {\left| {{{\tilde x}_{k - 1}}} \right.} \right\} - {{({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})}^T} \times }\\ &{\left[ {{H_k}} \right.({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})H_k^T}+\\ &{ {{\left. {{R_k}} \right]}^{ - 1}} \times \left({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}}\right)} \end{split}$$ (54)

      式(54)的最后一项, 根据引理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{ x}}_{k - 1}}).}\\[-10pt] \end{split}$$ (55)

      式(55)分别左右乘${\tilde x_{k - 1}}$

      $$\begin{split} &\tilde{ x}_{k - 1}^T\hat P_{k - 1}^{ - 1}{{\tilde x}_{k - 1}} > {({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})^T}[({H_k}{\beta _k}{F_k}){{\hat P}_{k - 1}} \times\\ & {({H_k}{\beta _k}{F_k})^T} + {H_k}{\Xi _k}H_k^T + {R_k}{]^{ - 1}}({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}}).\\[-10pt] \end{split}$$ (56)

      选择

      $$\begin{split} {\lambda _k} =& \Bigl\{ {({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})^{\rm T}}\Bigl[({H_k}{\beta _k}{F_k}){{\hat P}_{k - 1}} \times \\ &{({H_k}{\beta _k}{F_k})^{ T}} \Bigl]+ {H_k}{Q _k}H_k^{\rm T} + {R_k}{]^{ - 1}}\times\\ & ({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})\Bigl\}/(\tilde{ x}_{k - 1}^{\rm{T}}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}}) \end{split}$$ (57)

      由式(56)和式(57)可知, ${\lambda _k} < 1$,

      $$\begin{split} &{\lambda _k} \geqslant {p_{\min }}{({h_{\min }}{\beta _{\min }}{f_{\min }})^2} \\ & \big[{p_{\max }}{({h_{\max }}{\beta _{\max }}{f_{\max }})^2} + \\ &{{\hat q}_{\max }}h_{\max }^2 + {r_{\max }}{\big]^{ - 1}} = {\lambda _{\min }} > 0 \\[-10pt] \end{split} $$ (58)

      由式(57)得到

      $$\begin{split} &- {({H_k}{\beta _k}{F_k})^T}[({H_k}{\beta _k}{F_k}){{\hat P}_{k - 1}}{({H_k}{\beta _k}{F_k})^T}+\\ & {H_k}{{\hat Q}_k}H_k^T + {R_k}{]^{ - 1}}({H_k}{\beta _k}{F_k}{{\tilde{ x}}_{k - 1}})\leqslant -\\ & {\lambda _{\min }}{V_{k - 1}}({{\tilde{ x}}_{k - 1}}) \end{split}$$ (59)

      考虑式(54)中的其他项

      $$\begin{split} &{{\mu _k} = E\left\{ {\omega _k^T} \right.\left[ {{{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})}^{ - 1}}} \right.}-\\ & \quad \; \quad{ {H_k}^T\left( {{H_k}} \right.({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})H_k^T + }\\ &\quad \; \quad{{{\left. {{R_k}} \right)}^{ - 1}} \times \left. {{H_k}} \right]{\omega _k} + }\\ &\quad \; \quad{\upsilon _k^TR_k^{ - 1}{H_k}{{\hat P}_k}H_k^TR_k^{ - 1}{v_k}\left. {\left| {{{\tilde x}_{k - 1}}} \right.} \right\}} \end{split}$$ (60)

      考虑式(60)两边都是标量, 取迹不会改变等式, 所以根据

      $$tr(AB) = tr(BA)$$ (61)

      得到

      $$\begin{split} &{{\mu _k} = tr\left\{ {\left[ {{{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})}^{ - 1}}} \right.} \right. - {H_k}^T \times }\\ &\quad\quad\;{\left( {{H_k}} \right.({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})H_k^T + }\\ &\quad\quad\;{{{\left. {{R_k}} \right)}^{ - 1}}\left. {{H_k}} \right]\left. {{Q_k}} \right\} + tr[R_k^{ - 1}{H_k}{{\hat P}_k}H_k^T]\left| {{{\tilde x}_{k - 1}}} \right.\} } \end{split}$$ (62)

      根据引理2

      $$\begin{split} &{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})^{ - 1}} - \\ &{H_k}^T[{H_k}({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + \\ &{{\hat Q}_k})H_k^T + {R_k}{)^{ - 1}}{H_k}] > 0 \end{split}$$ (63)

      由式(63)可知${\mu _k} > 0$, 由式(42)和式(43)得

      $$\begin{split} {\mu _k} \leqslant & tr[{({\beta _k}{F_k}{{\hat P}_{k - 1}}F_k^T{\beta _k} + {{\hat Q}_k})^{ - 1}}{{\hat Q}_k}]+\\ & tr(R_k^{ - 1}{H_k}{{\hat P}_k}H_k^T)\leqslant\\ & \dfrac{{{q_{\max }}}}{{{{\hat q}_{\min }}}} \cdot L + \dfrac{{{h_{{{\max }^2}}}{p_{\max }}}}{{{r_{\min }}}} \cdot M = {\mu _{\max }} \end{split}$$ (64)

      因此根据引理1, 得到

      $$E[{V_k}({x_k})] - {V_{k - 1}}({\tilde x_{k - 1}}) \leqslant {\mu _{\max }} - {\lambda _{\min }}{V_{k - 1}}({\tilde x_{k - 1}})$$ (65)

      保证是${\tilde x_k}$有界的,

      $$\begin{split} E[{\left\| {{{\tilde{ x}}_k}} \right\|^2}] \leqslant& \dfrac{{{p_{\max }}}}{{{p_{\min }}}}E[{\left\| {{{\tilde{ 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{ x}}_k},k \geqslant 1} \right\}$是随机过程, $E({\left\| {{{\tilde{ x}}_k}} \right\|^2})$是该随机过程的均值函数, 随机有界指该随机过程有界. 即, 对于任意时刻$k$, 对于服从高斯分布的噪声$\omega $, 定理确保了误差${\tilde{ x}_k}$的2-范数平方的期望值是有界的, 即$E({\left\| {{{\tilde{ x}}_k}} \right\|^2}) \leqslant \alpha $上限$\alpha $的值由定理确定.

    • 该部分分析了文献[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 (\dfrac{\alpha }{{{p_t}}})$$
      $$\begin{split} & \Gamma = \inf \left\{ {t:{g_t} \geqslant 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} \leqslant 0$. 当${y_t}$是异常值时, ${g_t}$的值会随时间累积直到超过检测阈值$h$($h$的确定参考文献[6]), 此时发出警报. 该算法通过引入$h$, 避免了欧式检测算法那种仅考虑误检率的阈值设计方法, 在误检率和误检周期之间做到了平衡, 既降低了噪声的影响, 又可以尽快检测出攻击.

      但是, 该算法实质上是依赖与噪声的高斯分布特性, 即, 实质上还是${\chi ^2}$检测, 所以在遇到本文所研究的攻击时, 攻击前后的残差基本保持不变, 该算法是不可用的.

      文献[10]中, 针对具有隐蔽性的复数域攻击, 提出一种检测思想, 滤波器在不受攻击下工作和滤波器在错误数据注入攻击下工作会产生两个服从高斯分布的残差序列, 通过计算上述两个服从高斯分布的残差序列的巴氏相似性系数$DB$来判断攻击是否发生.

      当物理系统的基本动力学用线性系统和高斯统计特性的噪声进行建模时, 在正常操作期间, 残差序列${v^o}$遵循零均值高斯分布${v^o} = N({r^o},{S^o})$, 攻击发生后残差序列$v_k^a = 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} = N({m_1},{P_1})$, 攻击发生后残差序列$v_k^a = N({m_2},{P_2})$, 实数域的巴氏距离公式

      $$\begin{split} DB(v_k^o,v_k^a) = &\dfrac{1}{8}[{({m_1} - {m_2})^T}{P^{ - 1}}({m_1} - {m_2}) +\\ & \dfrac{1}{2}\ln \dfrac{{\det P}}{{\sqrt {\det {P_1}\det {P_2}} }} \\ \end{split} $$
      $$P = \frac{{{P_1} + {P_1}}}{2}$$

      Bhattacharyya距离无法辨别虚假数据注入攻击前后两个残差序列的相似性, 因为虚假数据注入攻击前后残差保持不变. 除此之外, 该检测法的前提是线性系统, 因为线性运算虽然会改变系统中随机变量的均值和方差但并不会改变随机变量服从高斯分布的特性, 对于本文所研究的非线性系统, 这种方法是不可用的.

      对文献[11], 目前存在的检测方法大多是在加性噪声的基础上提出的, 即假设过程方程和量测方程相对于噪声是线性的,

      $${{x}_k} = f({{x}_k},{u_k},{t_k}) + {\omega _k}$$
      $${{y}_k} = h({{x}_k},{t_k}) + {v_k}.$$

      文献[11]中欧式检测法的阈值和误检率正是基于噪声项的附加性得到的, 因为噪声并不参与到非线性变换中, 所以噪声的特性就是假设的正态分布的特性, 使用正态分布的$3\sigma $原则设计阈值. 只有在这种前提下, 欧氏距离检测算法才是正确的. 实际上, 过程和量测方程相对于噪声也可能是非线性的,

      $${{x}_k} = f({{x}_k},{u_k},{t_k},{\omega _k})$$
      $${y_k} = h({{x}_k},{t_k},{v_k}).$$

      在这种情况下, 欧氏距离检测法的阈值设计思想是行不通的, 因为噪声不再作为附加项, 而是直接参与到非线性变换中. 此时, 噪声不再服从正态分布, 统计特性未知, 阈值无法设计.

      综上, 本文提出一种检测思想: 在攻击未发生时, 利用本文提出的非线性滤波方法, 根据总线$i$上传感器在时间范围$[{t_0},{t_k}]$收集到的测量数据获得系统状态估计值, 进而获得后验状态误差值${\tilde{ x}_k}$, 此时我们并不知道状态误差${\tilde{ x}_k}$这个随机变量服从什么分布, 而且也不用对该随机变量的统计特性作出假设. 我们将时间区间$[{t_0},{t_k}]$分成$l$个小区间, 对每一个小区间包含的${\tilde x_k}$使用中心极限定理[18]构造随机变量${{S}_i}\tilde x = \frac{1}{T}\sum\nolimits_{k = 1}^T {{{\tilde{ x}}_k}}$, 得到$l$${{S}_i}\tilde x,i \in l$, 这$l$个随机变量是服从标准正态分布的, 当攻击发生时, ${S_i}\tilde x$会偏离标准整体分布, 从而达到检测攻击的目的.

      本文所提出的检测方法从系统内部状态误差入手, 不受限于传感器参数的概率分布未知所带来的影响, 不受限于加性噪声和非加性噪声的情况, 可以成功完成攻击检测, 并从数学角度给出了误检率与阈值之间的关系. 具体计算细节见下文.

    • 式(24)定义了随机变量${\tilde{ x}_k}$, 现在定义当系统未遭受攻击时, 记变量为${\tilde{ x}_k}$; 当系统遭受攻击时, 记变量为$\tilde{ x}_k^a$. 下面分别计算这两个随机变量的两个统计特性, 均值和方差.

      当系统正常运行时,

      $$\begin{split} E({{\tilde{ x}}_k}) =& E({{x}_k} - {{\hat{ x}}_{k\left| {k - 1} \right.}} - {K_k}({{y}_k} - {H_k}{{\hat{ x}}_{k\left| {k - 1} \right.}}))\\ = &E({{x}_k} \!-\! {{\hat{ x}}_{k\left| {k - 1} \right.}} \!-\! {K_k}({H_k}{{\hat{ x}}_k} \!+\! {v_k}\! -\! {H_k}{{\hat{ x}}_{k\left| {k - 1} \right.}}))\\ = &E(({I_k} - {K_k}{H_k})({{x}_k} - {{\hat{ x}}_{k\left| {k - 1} \right.}}) - {K_k}E({v_k}))\\ = &({I_k} - {K_k}{H_k})E({{x}_k} - {{\hat{ x}}_{k\left| {k - 1} \right.}}) - {K_k}E({v_k}) \end{split}$$ (67)

      当初值设置为${\hat{ x}_0} = E({{x}_0})$时, 上式变成

      $$E({\tilde{ x}_k}) = 0$$ (68)

      下面计算方差

      $$\begin{split} &{{{\hat P}_{k\left| k \right.}}\! =\! E({{\tilde{ x}}_k}\tilde{ x}_k^T)}\\ &\quad\quad\;{ \!=\! E\left\{ {\left[ {({I_k} - {K_k}{H_k})} \right.} \right.({{x}_k} - {{\hat{ x}}_{k\left| {k - 1} \right.}}) \!-\! \left. {{K_k}E({v_k})} \right]\! \times }\\ &\quad\quad\quad{{{[({I_k} - {K_k}{H_k})({{x}_k} - {{\hat{ x}}_{k\left| {k - 1} \right.}}) - {K_k}E({v_k})]}^T}}\\ &\quad\quad\;{\! =\! ({I_k} - {K_k}{H_k}){{\hat P}_{k\left| {k - 1} \right.}}{{({I_k} - {K_k}{H_k})}^T} + {K_k}{R_k}K_k^T} \end{split}$$ (69)

      当系统遭受攻击时,

      $$\begin{split} E(\tilde{ x}_k^a) &= E({x}_k^a \!-\! \hat{ x}_{k\left| {k - 1} \right.}^a \!-\! {K_k}({y}_k^a\! +\! {a_k} \!-\! {H_k}\hat{ x}_{k\left| {k - 1} \right.}^a))\\ &= E({x}_k^a - \hat{ x}_{k\left| {k - 1} \right.}^a - {K_k}\\ &\quad({H_k}{x}_k^a + {v_k} + {a_k} - {H_k}\hat{ x}_{k\left| {k - 1} \right.}^a))\\ &= E(({I_k} - {K_k}{H_k})({x}_k^a - \hat{ x}_{k\left| {k - 1} \right.}^a) - {K_k}E({v_k})) \end{split}$$ (70)

      下面计算方差

      $$\begin{split} \hat P_{k\left| k \right.}^a =& E[(\tilde{ x}_k^a){(\tilde{ x}_k^a)^T}]\\ =& E\{ [({I_k} \!- \!{K_k}{H_k})({x}_k^a \!-\! \hat{ x}_{k\left| {k \!-\! 1} \right.}^a) \!-\! {K_k}E({v_k} \!+\! {a_k})] \!\times\! \\ &{[({I_k} - {K_k}{H_k})({x}_k^a - \hat{ x}_{k\left| {k - 1} \right.}^a) - {K_k}E({v_k} + {a_k})]^T}\\ =& ({I_k} - {K_k}{H_k})\hat P_{k\left| {k - 1} \right.}^a{({I_k} - {K_k}{H_k})^T} + {K_k}{R_k}K_k^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{ x}}_k}}$, 根据中心极限定理, 可知

      $$S\tilde x \sim N\left(\frac{1}{T}\sum\limits_{k = 1}^T {E(\tilde{ 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{ 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 {E(\tilde{ 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 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{array}{*{20}{c}} {{H_0}:S \in [{T_l},{T_h}]} \\ {{H_1}:S \notin [{T_l},{T_h}]} \end{array}} \right.$$ (72)

      ${H_0}$代表系统没有遭受攻击, ${H_1}$代表系统遭受攻击.

      误检率${P_F}$的定义是攻击没有发生却宣布攻击发生的概率, 即$S\tilde x \sim N(m,\upsilon )$, 同时$\left| {S\tilde x - m} \right| > k\upsilon $.

      $$\begin{split} {P_F}(k) &= 1 - \int\limits_{m - kv}^{m + kv} {\frac{{\exp ( - {{(u - m)}^2}/(2{v^2}))}}{{v\sqrt {2\pi } }}} du \\ & = 1 - \Phi (k) + \Phi ( - k) \\ \end{split} $$
    • 本节验证本文设计的自适应平方根无迹卡尔曼滤波算法估计状态及所提出检测算法检测攻击有效性, 并与欧几里得检测方法作对比. 仿真参数和初始值如下, $m = [8.9,{{ }}8.8,{{ }}8.5]$, $d = [3.1,{{ }}3.4,{{ }}3.7]$ $u = [6.3, {{ }}1.6,{{ }}8.5]^T$, $B = [{0_{3 \times 3}};diag(1/m)]$.

      为了最小化由噪声引起的误报率, 欧几里得检测法的阈值设计采用$3\sigma $准则($\sigma $是噪声信号的标准差), 那么得到先验阈值$f = 3\sigma = 0.3$.

    • 对于模型(6), 未遭受攻击时的状态估计, 从图3可以看出, ASRUKF滤波器的滤波性能很好, 滤波精度很高.

      图  3  ASRUKF下的状态估计

      Figure 3.  State estimation in ASRUKF

    • 第一种检测方法, 欧几里得检测法文献[11]. 由图4可知, 在$k = 30$时刻对系统注入攻击, 在$k = 80$时刻成功检测到攻击, 欧几里得检测法可以成功检测出虚假数据注入攻击.

      图  4  两种检测方法针对隐蔽假数据攻击

      Figure 4.  Two detection methods for covert false data attack

      第二种检测方法, 巴氏相似性系数检测法文献[10].

      巴氏系数范围为0到1之间. 越接近0, 证明两序列越相似. 由图5可知, 巴氏相似性系数无法辨别虚假数据注入攻击前后两个残差序列的相似性, 不可以检测出虚假数据注入攻击.

      图  5  巴氏相似性系数

      Figure 5.  Bhattacharyya coefficient

      第三种检测方法是本文所提的检测方法. 首先, 使用分位数-分位数(Q-Q)图来确定前面根据中心极限定理构造的$S\tilde x$分布可以近似为高斯分布. 因为根据中心极限定理, 随机变量的数量越多, 它们的和也越近似服从高斯分布. 然而, 考虑实际情况, 取无限多的独立随机变量是不符合实际的, 所以在近似高斯分布的同时, 要尽量选取较少的独立随机变量. 要利用Q-Q图鉴别样本数据是否近似于高斯分布, 只需看Q-Q图上的点是否近似地在一条直线附近.

      本文的检测算法选取了$T = 30$个独立随机变量参与运算, 由图6可知, 数据的分布非常接近直线, 这证明$S\tilde x$的分布可以视为高斯分布.

      图  6  $S\tilde x$的Q-Q图

      Figure 6.  quantile-quantile plot of $S\tilde x$

      图7可以看出, 本文的检测算法对于隐蔽假数据攻击的有效性, 在$k = 30$处注入攻击, 选取30个迭代步骤下的随机变量构造出$S\tilde x$, 所以在$k = 60$处超过了阈值, 实现了攻击检测.

      图  7  本文提出的攻击检测方法

      Figure 7.  Attack detection proposed in this paper

      图8中给出误检率随着$k$变化的曲线.

      图  8  误检率${P_F}$

      Figure 8.  False alarm rate ${P_F}$

      图8可知, 一个算法的误检率过高, 是因为阈值设置的太低, 导致攻击之外的因素也会使得被检测的量超过阈值. 随着$k$变大, 阈值变高, 误检率自然就下降了.

      图8中, 当$k = {{3}}$时, 误检率${P_F} = 0.0027$, 除此之外, 由图4可知, 欧式检测法检测到攻击所用时间是$\Delta k = 50$; 由图8可知, 检测到攻击所用时间是$\Delta k = 30$, 检测用时更短.

      综上, 对比三种检测方法, 本文提出的检测方法可以成功检测攻击不受限于传感器参数统计特性未知所带来的影响, 不受限于加性噪声和非加性噪声的情况, 对于实际系统的适用性更强.

    • 本文研究了智能电网系统中虚假数据注入攻击的检测问题. 针对非线性系统, 噪声统计特性未知的情况, 本文使用自适应平方根无迹卡尔曼滤波算法对系统内部状态和噪声作出估计. 针对传感器参数统计特性未知和非加性噪声的情况, 利用中心极限定理构造出符合正态分布的随机变量, 基于该随机变量提出了一种攻击检测方法并从评价指标的角度对算法作出分析, 该算法对于实际系统的适用性更强.

WeChat 关注分享

返回顶部

目录

    /

    返回文章
    返回