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

本文链接：https://blog.youkuaiyun.com/F_Matlab/article/details/149219142

👨‍🎓个人主页

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

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

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

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

目录

💥1 概述

一、KF状态估计在电力系统中的应用基础

二、虚假数据注入攻击（FDIA）的原理与实施

（一）攻击基础模型

（二）针对KF的动态攻击特性

三、FDIA对KF状态估计的影响机制

四、防御机制研究进展

（一）检测技术

（二）主动防御

五、典型研究案例

六、未来挑战与方向

结论

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

💥1 概述

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

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

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

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

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

一、KF状态估计在电力系统中的应用基础

卡尔曼滤波（KF）通过动态融合系统模型与实时量测数据，实现电力系统状态的最优估计。其核心包括预测-更新迭代过程：

数学基础
非线性系统扩展
- EKF：通过泰勒展开线性化非线性系统，但忽略高阶项导致截断误差。
- CKF：基于球面径向准则，避免Jacobian矩阵计算，在IEEE 30/246节点系统中验证了更高精度。
- UKF：利用无迹变换，但参数选择影响稳定性。

应用场景

领域	具体应用
负荷预测	多节点系统短期/长期负荷需求预测，支持调度决策
电能质量监控	谐波失真、不平衡度分析，THD计算
发电机动态状态估计	功角、转速、电磁功率实时追踪（如新英格兰68节点系统）
电压稳定性评估	实时电网状态估计与崩溃预警

二、虚假数据注入攻击（FDIA）的原理与实施

（一）攻击基础模型

攻击定义：
攻击路径
- 量测篡改（Path 1）：入侵RTU/PMU通信链路修改上传数据
- 控制信号篡改（Path 2）：伪造控制中心指令
- GPS时间同步攻击：欺骗PMU时钟同步信号

（二）针对KF的动态攻击特性

KF脆弱性根源：依赖历史状态的递推特性使攻击具有累积效应。单次小幅度攻击在迭代中被放大，导致状态估计持续偏移。
案例验证：
在IEEE 14节点系统中，FDIA使线路1-2有功功率升至131%、无功功率升至364%，残

三、FDIA对KF状态估计的影响机制

影响维度	具体表现
状态估计偏差	发电机功角估计误差超10%，电压幅值偏移0.1346 p.u.，相角偏移7.2496°
系统控制失效	基于错误状态的调度指令引发线路过载、频率失稳
经济与安全风险	乌克兰电网事件表明FDIA可导致大范围停电，经济损失超$10亿

四、防御机制研究进展

（一）检测技术

残差-机器学习融合
- LSTM/GNN预测模型：利用状态预测误差检测异常（IEEE 118节点验证精度达92%）
自适应滤波改进
- AUKF算法：通过Holt指数平滑法增强迟滞性，WLS与AUKF估计值偏差超阈值即报警。
- 融合CKF：局部估计量替代全局估计，攻击定位误差<3%。

（二）主动防御

关键量测保护
- 最优保护问题转化为Steiner树问题，保护15%关键节点可使攻击成本提升300%。
多源数据融合
- PMU与SCADA数据交叉验证，降低单点篡改风险。
动态拓扑变换
- 随机调整网络拓扑结构（如开关重组），破坏攻击者对 HH 矩阵的认知。

五、典型研究案例

发电机动态FDIA（李扬等，2019）
- 攻击方法：基于量测残差不变原理设计 aa，在IEEE 9节点系统中成功规避RCKF检测。
KF-Matlab仿真平台（研学社，2023）
- 可视化攻击影响：状态估计误差较正常情况扩大5倍（见下图）。
  ![] 注：仿真结果展示攻击后状态跟踪偏离真实值
部分信息攻击（2014）
- 仅需30%拓扑信息即可构造有效 aa，降低攻击门槛。

六、未来挑战与方向

挑战
- 高维非线性系统攻击向量快速生成（如GAN强化攻击）
- 分布式能源接入扩大攻击面
研究方向
- 区块链量测数据存证（防篡改）
- 联邦学习框架下的协同防御

结论

KF状态估计的递推特性使其成为FDIA的高价值目标。防御需结合改进估计算法（如AUKF/CKF）、多源校验及动态防护策略，同时探索人工智能与区块链的前沿应用。当前研究已验证机器学习检测与拓扑重构的有效性，但实时性与大规模部署仍需突破。

📚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