NRD安全估计算法的仿真过程及思考

2021/9/21 17:40:36

编程Tag： 算法检测攻击 err 噪声仿真量测 NRD

本文主要是介绍NRD安全估计算法的仿真过程及思考，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

文章目录

算法思考
- 算法原理
仿真背景
- 模型设置
- 攻击设置（与RWD中完全相同）
代码
效果及分析
- 隐匿攻击
- - 分析
  - 噪声设置
  - 追踪效果
  - 估计误差
  - 量测值记录
  - 攻击检测
  - 检测率
- 非隐匿攻击
- - 分析
  - 噪声设置
  - 追踪效果
  - 估计误差
  - 量测值记录
  - 攻击检测
  - 检测率
- 极端非隐匿攻击
- - 分析
  - 噪声设置
  - 追踪效果
  - 估计误差
  - 量测值记录
  - 攻击检测
  - 检测率

算法思考

算法原理

本质上说，NRD其实就是一维的（标量的）卡方检测+标准Kalman滤波+恢复算法。因此其攻击检测原理与标准卡方检测一致，都是考量两种分布的一致性情况，一致性越差，结果参量越大。于是就可以通过这种方法检测出攻击下的分布和期望的分布的不同，进而检测出攻击的发生。
因此可以这样理解：
受到非隐匿攻击，量测值显而易见的会发生分布的变化，这可以被检测出来，无需多言。
受到了（零均值）隐匿攻击后，基于卡方检测的攻击检测依然是有效的，因为量测值同样可能会发生期望值的偏移，进而以一定概率引发分布的变化，使得攻击得以检测出来，从这个角度讲，NRD对隐匿攻击的检测的思想与RWD相同。

仿真背景

模型设置

与RWD的仿真类似，为了简单，使用一个不相关双通道两传感器的模型，两个传感器每时刻独立地观测两个状态量，并进行融合。
使用背景同样为智能电网，电网中的参数变化缓慢，一般可以认为是稳定系统，因此设定状态转移矩阵A为单位矩阵，系统参数波动仅由噪声引起。

攻击设置（与RWD中完全相同）

我们设置第2个传感器受到攻击，第1个传感器不受攻击。

攻击发生在30 ~ 50 和 70 ~ 90时刻。

同时为了考量虚警和漏警，我们进行多次重复试验。

依照文献设置，分为三种情况：
1、隐匿攻击
设置攻击功率（标准差）=量测噪声功率（标准差）=1.1
2、非隐匿攻击
设置攻击功率（标准差）=2.6
量测噪声功率（标准差）=1.1
3、极端非隐匿攻击
设置攻击功率（标准差）=5
量测噪声功率（标准差）=0.1

代码

%Project: 不相关&双通道&二维&新型残差检测器
%Author: Jace
%Data: 2021/09/21
%--------------------准备---------------------
close all;
clear all;
clc;
%------------------初始化参数---------------------

N=100;%设定采样点数，即持续时长

%攻击检测相关参数
H=3;%检测阈值

gk1=zeros(1,N);%检测量
gk2=zeros(1,N);%检测量

%噪声相关参数
P0=0.01;%初始状态噪声协方差
Q=[0.01,0.0;0.0,0.01];%设定系统噪声
R1=0.1;%设定局部观测1噪声
R2=0.1;%设定局部观测2噪声
w=sqrt(Q)*randn(2,N);
v1=sqrt(R1)*randn(1,N);
v2=sqrt(R2)*randn(1,N);

%攻击相关参数
attack1=0;%攻击1是否开启
attack2=1;%攻击2是否开启

atk1=sqrt(R1)*randn(1,N);%隐匿攻击1
atk2=sqrt(R2)*randn(1,N);%隐匿攻击2，通过设置此处来变更攻击类型

%系统模型参数
A=[1,0;0,1];%状态转移矩阵
H1=[1,0];%局部量测1量测矩阵
H2=[0,1];%局部量测2量测矩阵

%----------------初始化分配空间-------------------------
%真实状态值初始化
x=zeros(2,N);%物体真实状态值（分配空间），2*N维矩阵
x(:,1)=[10+P0*randn(1);20+P0*randn(1)];%物体初始真实状态值

%误差协方差初始化
p1=[1,0;0,1];%局部滤波1误差协方差初始值
p2=[1,0;0,1];%局部滤波2误差协方差初始值
p=[1,0;0,1];%全局滤波误差协方差初始值

%两个传感器量测值初始化
z1=zeros(1,N);%量测值1
z1(1)=H1*x(:,1);%观测1真实值初始值,第一列的列向量取出第一行，成为标量

z2=zeros(1,N);%量测值2
z2(1)=H2*x(:,1);%观测2真实值初始值,第一列的列向量取出第二行，成为标量

%各滤波器估计值初始化
xkf1=zeros(2,N);%局部估计状态1
xkf1(:,1)=x(:,1);%局部估计状态1初始化为第一列的列向量

xkf1NRD=zeros(2,N);%局部估计状态1备份
xkf1NRD(:,1)=xkf1(:,1);

xkf2=zeros(2,N);%局部估计状态2
xkf2(:,1)=x(:,1);%局部估计状态2初始化为第一列的列向量

xkf2NRD=zeros(2,N);%局部估计状态2备份
xkf2NRD(:,1)=xkf2(:,1);

xkfNRD=zeros(2,N);%全局估计状态
xkfNRD(:,1)=x(:,1);%全局估计状态初始化为第一列的列向量

I=eye(2);%2*2单位矩阵

Alarm=zeros(2,N);
%----------------NRD------------------------
for k=2:N
    %系统模型
    x(:,k)=A*x(:,k-1)+w(k);
    %量测模型，标量
    z1(k)=H1*x(:,k)+v1(k);
    z2(k)=H2*x(:,k)+v2(k);
    
    %注入攻击
    if attack1==1
        if k>=30&&k<=50||k>=70&&k<=90
            z1(k)=z1(k)+atk1(k);
        end
    end
    
    if attack2==1
        if k>=30&&k<=50||k>=70&&k<=90
            z2(k)=z2(k)+atk2(k);
        end
    end
          
    %----------------标准Kalman过程------------------------
    %估计误差协方差预测更新
    p_pre1=A*p1*A'+Q;
    p_pre2=A*p2*A'+Q;
    %增益矩阵预测更新
    K1=p_pre1*H1'/(H1*p_pre1*H1'+R1);
    K2=p_pre2*H2'/(H2*p_pre2*H2'+R2);
    %状态值预测更新
    x_pre1=A*xkf1NRD(:,k-1);
    x_pre2=A*xkf2NRD(:,k-1);
    %备份先验估计值，用于存在攻击时的状态恢复
    xbup1=x_pre1;
    xbup2=x_pre2;
    %状态值量测更新
    xkf1(:,k)=x_pre1+K1*(z1(k)-H1*x_pre1);
    xkf2(:,k)=x_pre2+K2*(z2(k)-H2*x_pre2);
    %估计误差协方差量测更新
    p1=(I-K1*H1)*p_pre1;
    p2=(I-K2*H2)*p_pre2;
    %局部滤波的估计值记录
    xkf1NRD(:,k)=xkf1(:,k);
    xkf2NRD(:,k)=xkf2(:,k);
    
    %----------------攻击检测NRD算法部分------------------------
    Sigrk1=H1*p1*H1'+R1;
    Sigrk2=H2*p2*H2'+R2;
    gk1(k)=(z1(k)-H1*xkf1(:,k))^2/Sigrk1;
    gk2(k)=(z2(k)-H2*xkf2(:,k))^2/Sigrk2;
    
    %----------------利用备份状态值重新估计,NRD恢复过程------------------------
    if gk1(k)>H
        Alarm(1,k)=1;
        z1back=H1*xkf1NRD(:,k-1);%量测值回退到上一时刻的估计值
        %局部滤波1的重新估计（利用备份）
        xkf1NRD(:,k)=xbup1+K1*(z1back-H1*xbup1);
    end  
    if gk2(k)>H
        Alarm(2,k)=1;
        z2back=H2*xkf2NRD(:,k-1);%量测值回退到上一时刻的估计值
        %局部滤波2的重新估计（利用备份）
        xkf2NRD(:,k)=xbup2+K2*(z2back-H2*xbup2);
    end  
    
    %----------------融合部分------------------------
    %两个局部滤波不相关条件下的简单全局融合
    p=inv(inv(p1)+inv(p2));%融合后的估计误差协方差
    xkfNRD(:,k)=p*(inv(p1)*xkf1NRD(:,k)+inv(p2)*xkf2NRD(:,k));%融合后的状态估计值
end

%--------------------------误差计算--------------------------------
%初始化
messure_err_x1=zeros(1,N);
messure_err_x2=zeros(1,N);
kalman_err_x1=zeros(1,N);
NRD_err_x1=zeros(1,N);
kalman_err_x2=zeros(1,N);
NRD_err_x2=zeros(1,N);
for k=1:N
    %量测误差
    messure_err_x1(k)=abs(z1(k)-x(1,k));
    messure_err_x2(k)=abs(z2(k)-x(2,k));
    %Kalman估计偏差
    kalman_err_x1(k)=abs(xkf1(1,k)-x(1,k));
    kalman_err_x2(k)=abs(xkf2(2,k)-x(2,k));
    %NRD估计偏差
    NRD_err_x1(k)=abs(xkfNRD(1,k)-x(1,k));
    NRD_err_x2(k)=abs(xkfNRD(2,k)-x(2,k));
end

%噪声图
figure;
subplot(2,2,1)
plot(w(1,:));xlabel('采样时间');ylabel('噪声');
title('第1状态值过程噪声');
subplot(2,2,2)
plot(w(2,:));xlabel('采样时间');ylabel('噪声');
title('第2状态值过程噪声');
subplot(2,2,3)
plot(v1);xlabel('采样时间');ylabel('噪声');
title('第1状态值第1传感器量测噪声');
subplot(2,2,4)
plot(v2);xlabel('采样时间');ylabel('噪声');
title('第2状态值第2传感器量测噪声');

%局部滤波器1
figure;
t=2:N;
plot(t,x(1,t),'-k.',t,z1(t),'-b.',t,xkfNRD(1,t),'-r.',t,xkf1(1,t),'-g.');
legend('第1状态值','第1量测值','第1状态量NRD估计值','第1状态量Kalman估计值');
xlabel('采样时间');ylabel('位置');
title('局部滤波器1跟踪状态');

%局部滤波器2
figure;
t=2:N;
plot(t,x(2,t),'-k.',t,z2(t),'-b.',t,xkfNRD(2,t),'-r.',t,xkf2(2,t),'-g.');
legend('第2状态值','第2量测值','第2状态量NRD估计值','第2状态量Kalman估计值');
xlabel('采样时间');ylabel('位置');
title('局部滤波器2跟踪状态');

%局部滤波器1误差
figure;
hold on,box on;
plot(messure_err_x1,'-b.');
plot(kalman_err_x1,'-g.');
plot(NRD_err_x1,'-r.');
legend('量测','Kalman估计','NRD估计');
xlabel('采样时间');ylabel('误差');
title('局部滤波器1误差');

%局部滤波器2误差
figure;
hold on,box on;
plot(messure_err_x2,'-b.');
plot(kalman_err_x2,'-g.');
plot(NRD_err_x2,'-r.');
legend('量测','Kalman估计','NRD估计');
xlabel('采样时间');ylabel('误差');
title('局部滤波器2误差');

%量测值记录
figure;
hold on,box on;
plot(z1);
plot(z2);
xlabel('采样时间');ylabel('数值');
legend('量测1','量测2');
title('量测值记录');

%检测记录
figure;
t=2:N;
subplot(2,1,1)
plot(t,gk1(t),'-b.',t,Alarm(1,t),'-r.');
legend('攻击检测衡量值1','是否存在攻击1');
title('检测攻击1');
subplot(2,1,2)
plot(t,gk2(t),'-b.',t,Alarm(2,t),'-r.');
legend('攻击检测衡量值2','是否存在攻击2');
title('检测攻击2');