针对KF状态估计的电力系统虚假数据注入攻击研究（Matlab代码实现）-CSDN博客

本文链接：https://blog.csdn.net/qq_46009436/article/details/146640329

👨‍🎓个人主页

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

一、引言

二、电力系统状态估计与KF

三、虚假数据注入攻击（FDIA）

四、针对KF状态估计的FDIA研究

五、结论与建议

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

💥1 概述

如今现代电力系统正在向智能化方向发展，大量的智能设备，如智能仪表和传感器，促进了电力系统在发电、变电、输电和配电模式方面的转变，使得智能电网成为一个典型的网络物理系统，即将物理电力传输系统和计算机网络相结合。在智能电网中，监督控制和数据采集系统（SCADA）实时收集外场设备通过网络发送来的数据，进行分析后向控制中心汇报收集到的信息，控制中心根据这些信息对电网的发电配电进行调整。

在享受智能电网带来便利的同时，由于使用了大量智能设备以及通过网络来收发数据，使得智能电网更容易遭到攻击者的攻击。最为显著的是通过有预谋地篡改智能设备的数据来达到攻击电力系统的目的，称之为虚假数据注入攻击（FDIA）。为了保障电力服务的平稳运行，智能电网的检测至关重要。

FDIA 的攻击目标主要包括测量单元、通信网络和控制设备。由于控制设备往往保护级别较高，而难以入侵，因此，FDIA 通常以前两种方式实施。在实施 FDIA 过程中，攻击者首要任务是入侵系统网络。例如，为了破坏测量单元（例如，远程终端单元和向量测量单元）的数据，攻击者往往利用加密和认证机制固有的漏洞来进行入侵并修改原始数据。在文献[17]中，研究人员揭示了通过欺骗向量测量单元全球定位系统的时间同步攻击。由于全球定位系统信号并无加密或授权机制，攻击者可以产生伪冒的全球定位系统信号，而接收器无法从原始数据中辨别出来。

攻击者在侵入系统网络后，可以获得测量修改权限。根据潮流规律，某一测量值的变化会引起相邻测量值的变化。当伪造一个数据（例如节点或传输线路的测量值）时，为了绕过坏数据检测，攻击者会考虑潮流定律，找到测量数据相应改变的最小空间。而将需要协同攻击的最小空间定义为该数据在执行 FDIA 时满足潮流规律和最优资源利用规律的最小相关空间[18]。

状态估计方法分为静态估计和动态估计。传统的电力系统状态估计方法采用的是加权最小二乘估计器等静态估计方法。这是建立在电力系统具有足够冗余且处于稳态状态的假设之上的[8]。然而，由于用电需求和发电量是实时变化的，因此真实的电力系统并不是在稳定状态下运行[19]。同时静态估计方法仅考虑当前时刻数据，没有关联先前的状态，因此其估计结果精度不佳。为了解决这些问题，动态状态估计器如卡尔曼滤波器被引入电力系统状态估计[20]。接下来讨论了不同的基于状态估计的 FDIA 检测方案，其分类如图 3.2 所示。（本文选择KF）

一、引言

现代电力系统正在向智能化方向发展，智能电网作为物理电力传输系统和计算机网络的结合体，其安全性和稳定性至关重要。然而，由于智能电网使用大量智能设备并通过网络收发数据，因此更容易受到攻击者的攻击。虚假数据注入攻击（FDIA）是一种通过有预谋地篡改智能设备的数据来达到攻击电力系统目的的攻击方式。本文旨在研究针对KF状态估计的电力系统虚假数据注入攻击，以提出有效的检测与防御策略。

二、电力系统状态估计与KF

状态估计方法：
- 状态估计方法分为静态估计和动态估计。
- 传统的电力系统状态估计方法采用加权最小二乘估计器等静态估计方法，但该方法仅考虑当前时刻数据，未关联先前的状态，因此估计结果精度不佳。
- 动态状态估计器如卡尔曼滤波器（KF）被引入电力系统状态估计，以克服静态估计方法的局限性。
KF基本原理：
- KF算法是一种递推最优估计理论，采用状态空间描述法对系统特性加以限定，通过递推形式实现状态预测过程。
- KF算法分为状态估计与更新两个环节，基于上一时刻的状态预测值以及系统状态转移方程，实现当前时刻的状态估计及协方差估计，并根据观测值进行状态更新。

三、虚假数据注入攻击（FDIA）

FDIA的定义与特点：
- FDIA是通过有预谋地篡改智能设备的数据来达到攻击电力系统目的的攻击方式。
- 攻击者通常利用系统漏洞入侵网络，获得测量修改权限，并伪造数据以绕过坏数据检测。
FDIA的攻击目标：
- FDIA的攻击目标主要包括测量单元、通信网络和控制设备。由于控制设备保护级别较高，难以入侵，因此FDIA通常以前两种方式实施。
FDIA的攻击过程：
- 攻击者首先入侵系统网络，获得测量修改权限。
- 根据潮流规律，伪造一个数据（如节点或传输线路的测量值）时，攻击者会考虑潮流定律，找到测量数据相应改变的最小空间，即满足潮流规律和最优资源利用规律的最小相关空间。

四、针对KF状态估计的FDIA研究

攻击模型：
- 建立针对KF状态估计的FDIA攻击模型，分析攻击者如何通过伪造数据来影响KF的状态估计结果。
仿真实验：
- 利用MATLAB等工具进行仿真实验，模拟FDIA对KF状态估计的影响。
- 通过对比真实状态与估计状态之间的差异，评估FDIA的攻击效果。
检测结果与分析：
- 分析仿真实验结果，揭示FDIA对KF状态估计的潜在威胁。
- 探讨如何改进KF算法或提出新的检测策略来抵御FDIA的攻击。

五、结论与建议

结论：
- 本文研究了针对KF状态估计的电力系统虚假数据注入攻击，通过仿真实验揭示了FDIA对KF状态估计的潜在威胁。
建议：
- 加强智能电网的安全防护，提高系统对FDIA的检测与防御能力。
- 深入研究KF算法及其改进方法，提高状态估计的准确性和鲁棒性。
- 加强与网络安全领域的合作与交流，共同应对智能电网面临的安全挑战。

综上所述，针对KF状态估计的电力系统虚假数据注入攻击研究对于保障智能电网的安全运行具有重要意义。未来应继续深化相关研究，提出更加有效的检测与防御策略。

📚2 运行结果

clf
subplot(311), plot(t,yt,'b',t,y_hat,'r--'),
xlabel('Number of Samples'), ylabel('Output')
title('Kalman Filter Response')
legend('True','Filtered')
subplot(312), plot(t,yt-y,'g',t,yt-y_hat,'r--'),
xlabel('Number of Samples'), ylabel('Error')
legend('True - measured','True - filtered')

subplot(313), plot(t, state,'g',t,[x1_hat, x2_hat,x3_hat],'r--'),
xlabel('Number of Samples'), ylabel('State')
legend('True - state','True - filtered')

....

figure;
plot(t, filtered, t, v);
ylim([-2 2]);
xlabel('time (in seconds)');
legA = legend('filtered','real');
set(legA,'FontSize',12)
title('signal vs. estimation');
grid on

[y_o,t_o,x_o] = lsim( sys, [u 0*w], t); % IDEAL behavior for comparison
figure(1)
for k = 1:length(sys.A)
subplot( round((length(sys.A)+1)/2 ),2,k)
plot(t,x_o(:,k),'k','linewidth',2); hold on;
plot(t,x_(:,k),'bo:','linewidth',2);
stairs(t,x(:,k),'r','linewidth',1.5);
xlabel('Time (s)'); ylabel('Amplitude'); title(['State x_' num2str(k) '(t)']); grid on;
legend('ideal','real','estimated'); hold off;
end
subplot( round((length(sys.A)+1)/2 ),2,k+1)
plot(t,y_o,'k','linewidth',2); hold on;
plot(t,z,'bo:','linewidth',2);
stairs(t,y,'r','linewidth',1.5);
xlabel('Time (s)'); ylabel('Amplitude'); title('Output y(t)'); grid on;
legend('ideal','measured','estimated');

%Generate process noise and sensor noise vectors using the same noise covariance values Q and R that you used to design the filter.
rng(10,'twister');
w1 = sqrt(Q)*randn(length(t),1);
v1 = sqrt(R)*randn(length(t),1);

%Finally, simulate the response using lsim.
out = lsim(SimModel,[u,w1,v1]);
%lsim generates the response at the outputs yt and ye to the inputs applied at w, v, and u. Extract the yt and ye channels and compute the measured response.

yt = out(:,1); % true response
ye = out(:,2); % filtered response
y = yt + v1; % measured response

%Compare the true response with the filtered response.
clf
subplot(211), plot(t,yt,'b',t,ye,'r--'),
xlabel('Number of Samples'), ylabel('Output')
title('Kalman Filter Response')
legend('True','Filtered')
subplot(212), plot(t,yt-y,'g',t,yt-ye,'r--'),
xlabel('Number of Samples'), ylabel('Error')
legend('True - measured','True - filtered')

%covariance of the error before filtering (measurement error covariance)
MeasErr = yt-yt;
MeasErrCov = sum(MeasErr.*MeasErr)/length(MeasErr) %=0
% and after filtering (estimation error covariance)
EstErr = yt-ye;
EstErrCov = sum(EstErr.*EstErr)/length(EstErr)
% so the Kalman filter reduces the error yt - y due to measurement noise

% Kalman filter
[y_hat, ~, ~, P_last, K_last, x_hat, K_steps, P_steps] = kalmanfilter(y, Q, R, StopTime);

e_real = x - x_hat; %state error

figure();
plot(t,y, t, y_hat);
ylim([-1.5,1.5]);
legA = legend('real','estimated');
title('true signal and filtered signal without attack');grid on;

部分代码：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% kalman filter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n = 1:length(t)-1
% prediction
x(:,k) = F*x(:,k-1) +G*u(k-1) +J*w(:,k-1); % 'a priori' predicted state
P = F*P*F' +J*Q*J'; % 'a priori' predicted covariance

% correction
K = P*H'*inv( H*P*H' +R ); % computes optimal gain
e(k) = z(k) -H*x(:,k); % measurement residual
x(:,k) = x(:,k) +K*e(k); % 'a posteriori' state estimative
P = (I -K*H)*P; % 'a posteriori' covariance
y(k) = H*x(:,k); % estimated output
k = k +1;

% update
K_m = [K_m K]; % stores current K
P_m = [P_m P]; % stores current P
end
x = transpose(x);
end